MANEUVER  DESIGN  EOR  EAST 
SATELLITE  CIRCUMNAVIGATION 

THESIS 


Stanley  D.  Straight,  Captain,  USAF 


AFIT/GA/ENY/04-M05 


DEPARTMENT  OF  THE  AIR  FORCE 
AIR  UNIVERSITY 

AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 


Wright-Patterson  Air  Force  Base,  Ohio 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official 
policy  or  position  of  the  United  States  Air  Force,  Department  of  Defense,  or  the  United 
States  Government. 


AFIT/GA/ENY/04-M05 


MANEUVER  DESIGN  FOR  FAST  SATELLITE  CIRCUMNAVIGATION 

THESIS 


Presented  to  the  Faculty 
Department  of  Aeronautics  and  Astronautics 
Graduate  School  of  Engineering  and  Management 
Air  Force  Institute  of  Technology 
Air  University 

Air  Education  and  Training  Command 
In  Partial  Fulfillment  of  the  Requirements  for  the 
Degree  of  Master  of  Science  in  Astronautical  Engineering 


Stanley  D.  Straight,  BS 
Captain,  USAF 

March  2004 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


AFIT/GA/ENY/04-M05 


MANEUVER  DESIGN  FOR  FAST  SATELLITE  CIRCUMNAVIGATION 


Stanley  D.  Straight,  BS 
Captain,  USAF 


Approved: 


_ //Signed// _  _ 

Steven  G.  Tragesser  (Chairman)  date 


_//Signed// _  _ 

William  E.  Wiesel  (Member)  date 


//Signed//  _ 

Robert  A.  Howard  (Member) 


date 


AFIT/GA/ENY/04-M05 


Abstract 

The  feasibility  of  satellite  operations  in  close  proximity  to  a  reference 
satellite  is  of  interest  for  both  civilian  and  military  applications.  One  such  operation  is 
circular  circumnavigation  in  a  time  period  less  than  the  orbital  period  of  the  reference 
satellite.  This  thesis  investigates  a  guidance  scheme  for  such  maneuvers  involving 
impulsive  burns  at  specific  points  within  a  specified  toroidal  region  centered  on  the 
circular- orbiting  reference  satellite.  Two  analytical  methods  for  determining  the 
magnitude  and  direction  of  the  impulses  are  demonstrated.  These  methods  are  then  used 
as  initial  estimates  in  an  optimization  scheme  to  produce  the  minimum  total  required 
impulse. 
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MANEUVER  DESIGN  EOR  EAST  SATELLITE  CIRCUMNAVIGATION 


I.  Introduction 

Operations  between  two  spacecraft  in  orbit  have  become  of  increasing  interest  to 
both  the  civilian  and  military  communities.  Satellite-to-satellite  operations  have  been 
demonstrated  since  the  beginning  of  human  space  endeavors.  The  most  common  relative 
satellite-to-satellite  (relative  motion)  operations  have  been  rendezvous  between  two 
cooperating  spacecraft,  but  other  proximity  operations  are  becoming  more  important.  In 
recent  times,  there  has  been  considerable  interest  in  orbiting  satellites  in  close  relative 
formations. 

Considerable  work  has  been  done  in  the  area  of  satellite  formation  flying:  the 
design  of  formations  (11;  13),  their  reconfiguration  and  maintenance  (5),  and  formation 
evolution  through  orbital  perturbations  (3;  12;  14;  15).  Portions  of  this  work  have 
focused  on  natural  motion  formations,  establishing  a  relative  position  and  velocity  with 
respect  to  a  reference  orbit,  and  allowing  the  natural  dynamics  to  produce  elliptical 
motion  in  the  relative  frame.  Reconfiguration  of  formations,  another  topic  of  study,  has 
focused  on  optimizing  propellant  expenditure  from  one  formation  to  another  without 
necessarily  focusing  on  the  shape  or  time  variation  of  the  flight  path.  (5;  13) 

Other  proximity  operations  maneuvers  are  becoming  more  important  in  planning 
for  such  missions  as  on-orbit  repair  and  refueling  as  well  as  potential  damage  inspection 
or  identification  of  Resident  Space  Objects  (RSO).  (7;  10)  Circumnavigating  a  chief 
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satellite  with  a  deputy  satellite  provides  the  ability  to  inspect  the  chief  from  a  variety  of 
viewpoints.  Often  these  viewpoints  are  required  to  be  at  a  constant  distance  from  the 
chief  and  therefore  requiring  circular  circumnavigations.  A  circumnavigation  is  defined 
as  the  deputy’s  flight  about  a  desired  circular  path  (nominal  path)  with  a  specified 
orientation  about  a  chief  spacecraft. 

Previous  Work 

The  minimum  propellant  required  for  a  circular  circumnavigation  is  the  natural 
motion  circular  formation  (11:7-8),  which  requires  the  initial  conditions  to  be  set  up  in  a 
very  specific  manner.  These  natural  motion  circumnavigations  all  have  a 
circumnavigation  period  (rotating  around  the  chief  satellite  through  360°)  equal  to  the 
period  of  the  chief.  This  period  is  on  the  order  of  90  minutes  for  Low  Earth  Orbiting 
satellites  with  an  altitude  around  400  km,  and  increases  as  the  altitude  increases. 

Circumnavigation  times  less  than  the  chiefs  orbital  period,  are  termed  ‘fast’. 
These  fast  circumnavigation  maneuvers  have  utility  in  the  operation  of  a  proposed 
‘inspector’  micro- satellite  (4:1).  In  order  to  determine  time  changing  phenomena  on  the 
chief  spacecraft,  these  maneuvers  need  to  be  accomplished  less  than  the  orbital  period  of 
the  chief  satellite  (4:1).  As  the  required  time  for  a  circumnavigation  decreases,  the  total 
impulse  required  increases  to  perform  the  circumnavigation.  Minimizing  the  total 
impulse  for  a  given  maneuver  allows  for  greater  operational  flexibility,  increased  sorties 
for  a  given  amount  of  propellant,  and  potentially  increase  the  lifetime  for  any  given 
satellite  performing  these  maneuvers. 
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The  theoretical  impulse  requirements  for  differing  circumnavigation  times  and 
orientations  of  a  nominal  circular  path  have  been  demonstrated  (4:  1-2).  Their  method 
assumed  continuous  control  to  produce  perfect  circular  motion  with  respect  to  the  chief. 
After  simplifications,  the  unperturbed  Hill’s  equations  are  derived  as  (4:2;  14:377) 

x  -  2ny  -  3n^x  — 

y  +  2nx^f^  (1) 

z  +  n^z=f, 

where  x,  y,  and  z  represent  the  position  as  a  function  of  time  relative  to  the  chief,  and  “fx, 
fy,  and  fz  are  the  propulsive  forces  per  unit  mass”  (4:2).  From  these  equations  the  total 
impulse,  required  to  follow  a  path  defined  using  Equation  1  is  represented  as  (4:2): 

T 

(2) 

0 

where  T  is  the  time  required  to  follow  the  total  path.  This  theoretical  impulse  is 
informative  in  understanding  the  total  forces  required  to  follow  a  specific  path,  but  often 
the  path  is  only  constrained  to  lie  within  a  certain  volume.  Exact  adherence  to  a  defined 
circumnavigation  path  is  not  always  necessary;  operationally,  it  is  conjectured  only  a 
very  small  percentage  of  the  total  flight  path  is  required  to  be  at  a  certain  position  and 
time  relative  to  the  chief.  The  rest  of  the  flight  path  is  then  constrained  to  be  distance 
away  from  nominal  path,  or  in  some  instances  a  minimum  distance  to  the  chief 
spacecraft/object. 
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Problem  Statement 


For  the  problem  investigated  here,  the  deputy  will  perform  the  circumnavigation 
through  a  set  of  discrete  impulses  in  a  specified  time  of  flight  (TOF).  The  deputy’s  flight 
path  must  perform  a  full  27U  angular  rotation  about  the  chief  without  doubling  back  on 
itself.  The  placement  and  number  of  each  bum  point  must  be  determined  as  well  as  the 
required  direction  and  magnitude  of  each  individual  impulse. 

The  deputy  is  also  constrained  to  stay  within  a  specified  distance  from  the  desired 
or  ‘nominal’  path  during  the  circumnavigation  maneuver.  This  constraint  allows  for 
operational  considerations  such  as  collision  avoidance  and  operational  requirements  for 
the  deputy’s  payload.  The  payload  is  postulated  to  potentially  be  a  remote  sensing 
detector  (visual,  infrared,  etc. . .)  where  the  distance  to  the  chief  can  become  an  important 
operational  factor. 

There  are  two  probable  cases  of  general  mles  on  constraining  the  placement  of  the 
burn  points.  First,  a  case  is  defined  by  requiring  all  the  bum  points  placed  on  the  nominal 
path.  This  is  called  the  ‘Special  Case’.  Next,  a  case  is  defined  by  allowing  the  burn 
points  to  be  placed  anywhere  within  the  constraint  volume,  called  the  ‘General  Case’. 

The  special  case  may  be  required  if  the  spacecraft  is  required  to  be  a  constant  distance 
from  the  chief.  For  instance,  there  may  be  a  plume  exclusion  zone  or  distance 
requirement. 

The  general  case  has  the  most  operational  utility,  because  remote  measurements 
of  the  chief  will  be  required  after  the  relatively  dynamic  behavior  of  the  spacecraft  settles 
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after  each  subsequent  burn  point.  Because  the  burn  points  do  not  have  to  be  on  the 
nominal  path,  it  allows  more  flexibility  in  choosing  where  bums  are  placed. 

Overview  of  Content 

A  method  involving  discrete  impulsive  burns  is  desired.  These  maneuvers  require 
significant  propellant  expenditure  to  achieve  the  required  circumnavigation  trajectories 
within  the  desired  TOP.  The  placement,  relative  to  the  chief  satellite,  and  the  timing  of 
these  discrete  impulses  (at  the  burn  points)  has  a  significant  impact  on  the  propellant 
required  for  a  given  maneuver. 

Hill’s  equations  are  used  as  the  primary  tool  for  modeling  the  dynamics  for  the 
required  maneuvers.  The  equations  are  used  to  calculate  the  total  impulse,  Avt,  required 
for  a  given  circumnavigation  maneuver.  From  these  equations,  the  position,  magnitude, 
and  direction  of  each  discrete  impulse,  Av,  is  determined.  The  magnitudes  are  summed 
to  determine  the  total  impulse  required.  The  required  impulse  is  directly  proportional  to 
the  amount  of  thrust  a  propulsion  system  must  provide,  and  thus  the  mass  of  propellant 
required  for  a  given  mission. 

A  simple  method  for  placement  of  the  bum  points  is  initially  developed,  and  used 
as  an  initial  guess  for  numerical  optimization.  The  optimization  routine  is  used  to 
investigate  the  lowest  minimum  propellant  required.  From  analyzing  the  optimization 
results,  an  analytical  algorithm  is  proposed  to  approximate  the  minimum  total  required 
impulse  for  a  circumnavigation  maneuver,  and  to  develop  a  more  robust  initial  guess  for 
the  numerical  optimization.  This  algorithm’s  performance  is  then  evaluated  for  several 
cases. 
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II  Methodology 


Assumptions 

The  chief  or  RSO  is  assumed  to  be  in  a  circular  orbit,  with  the  deputy  orbit  having 
a  very  small  eccentricity.  Additionally,  a  two  body,  point- mass  gravitational  model  is 
assumed  (i.e.  no  perturbations),  and  the  distance  between  the  deputy  and  the  chief  is 
much  less  than  the  radius  of  the  chiefs  orbit.  These  assumptions  allow  the  use  of  Hill’s 
equations  for  relative  orbital  motion. 


Hill’s  Equations. 


A  specific  form  of  Hill’s  equations  (16:83)  is  used  (also  known  as  the  complete 
Clohessy- Wiltshire  solution)  and  shown  as: 


5r(r) 

5v(0 


Sv(to). 


(3) 


The  matrix  (O)  (16:83)  are  defined  as: 


O(A0  = 


rr  n 


4-3cos(/iAO  0  0 

6(sm(nAt)-  nAt)  1  0 


0 


3nsin(nAt)  0 

6n(cos{nAt)-l)  0 
0  0 


—  sin(/iA0  —  (1 -cos(/iAr)) 

n  n 

2  4 

—  (cos(fzAO-l)  —sin(nAt)  -3nA: 


0  cos(/iA0 


0 

0 

-n  sin(nAt) 


0 

cos(nAt) 

-2sin(nAt) 

0 


0 

2sin(nAt) 

-3  +  4cos(/iA0 
0 


0 

0 

isin(/iA0 

n 

0 

0 

cos(/iA0 


(4) 


where  n  is  the  mean  motion  of  the  chief’s  orbit  and  At=t-to.  Equation  3  determines  the 
position,  5r,  and  velocity,  6v,  relative  to  the  chief  at  a  time.  At,  later  than  the  initial 
position  and  velocity. 
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Assumptions  in  Initial  Conditions. 


The  deputy’s  initial  position  and  velocity  relative  to  the  chief  are  assumed  to  be 
known.  However,  in  order  to  generalize  the  results,  zero  initial  relative  velocity  is  used 
for  all  calculations.  Zero  initial  velocity  ensures  no  component  direction  of  the  initial 
velocity  can  subtract  or  add  to  the  impulse  at  the  initial  point.  This  essentially  cancels  the 
effects  of  any  variations  in  the  initial  conditions  on  the  overall  optimization. 

Additionally,  the  deputy  is  required  to  not  exceed  a  given  distance  from  the 
nominal  path  during  any  part  of  the  maneuver.  The  distance  from  the  deputy  to  the 
nominal  path  is  defined  as  a  maximum  deviation,  pmax-  This  deviation,  Pmax,  is  measured 
as  the  magnitude  of  the  spatial  deviation  vector  of  the  flight  path  from  the  nominal  path. 
This  spatial  deviation  vector  is  thought  of  as  a  radius  from  the  nominal  path,  thus  pmax  is 
termed  the  maximum  deviation  radius.  The  actual  deviation  (of  the  flight  path  from  the 
nominal  path)  is  only  measured  in  a  spatial  sense.  It  does  not  take  into  account  when  and 
where  the  deputy  is  located  on  the  flight  path  with  reference  to  when  and  where  the 
deputy  is  to  be  nominally  located  along  the  nominal  path.  The  pmax  constraint  defines  a 
toroidal  constraint  surface  about  the  nominal  path. 

Instantaneous  Av  Assumption. 

Instantaneous  impulses,  Av,  which  occur  at  discrete  points  in  space,  are  assumed. 
This  assumption  is  less  valid  for  low  thrust  vehicles  or  for  extremely  fast 
circumnavigation  times  of  flight  when  impulses  may  require  a  significant  amount  of  time 
to  impart  a  change  in  velocity.  For  instance,  this  assumption  breaks  down  as  the 
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individual  maneuver  durations  become  a  significant  fraction  of  the  circumnavigation  time 


of  flight,  TOF. 

Nominal  Path 


The  orientation  and  size  of  the  circular  nominal  path  can  be  described  by  four 
parameters:  ro,  ?,  Ty,  and  Tz.  (4:3)  A  2-3  space  fixed  Ty,  Tz  rotation  sequence  of  a  circle 
of  radius,  ro,  in  the  y-z  plane  defines  the  nominal  circular  path.  The  angle  ?  defines  a 
spatial  degree  of  freedom  along  the  circular  path  with  the  initial  point  being  defined  byyo. 
Figure  1  illustrates  this  rotation.  The  values  of  Ty,  Tz,  and  ro  are  assumed  to  be  given 
quantities,  whereas  ?  is  a  variable  which  must  be  varied  to  determine  points  on  the  circle. 

These  four  parameters  are  defined  in  the  Local  Vertical,  Local  Horizontal  (LVLH) 
coordinate  system.  The  LVLH  coordinates  define  the  y  direction  in  the  same  direction  of 
the  chiefs  instantaneous  velocity  vector.  The  x  direction  is  defined  in  the  radial  (from 
the  center  of  the  Earth)  direction  to  the  chief,  and  consequently  the  z  direction  is 
orthogonal  to  x  and  y.  This  coordinate  system  is  equivalent  to  the  RSW  coordinates  used 
in  many  texts  (14:162-163). 


Figure  1  a)  Rotated  Nominal  Path  b)  Unrotated  Nominal  Path 
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Mathematically,  any  point  along  the  nominal  path  is  represented  by  a  phasing 


angle,  Yi,  and  by  rotating  an  initial  vector  of  length  ro  placed  along  the  y  direction  ([0;  ro; 
0]).  This  corresponds  to  the  initial  point  located  with  Jq  =  0.  Using  the  rotation  described 
above,  the  position  vector  as  a  function  of  y  is  (4:3) 


6r(Y,) 


Tq  •  cos(0j, )  •  sin(  0y )  •  sin(  Y,)  -  Tq  •  sin(  0^ )  •  cos(y,) 
Tq  •cos(0j)-cos(y,)  +  Tq  •  sin(Y,)-sin(  0j)-sin(  0^) 
ro-cos(0y)-sin(Y,) 


(5) 


where  Yo  =  0. 


State  Vector  Deflnition 

Unique  spatial  positions  where  discrete,  instantaneous  impulses  occur  are 
called  ‘burn  points’.  These  burn  points  are  required  to  perform  the  circumnavigation 
within  the  required  total  time.  Assigning  individual  time  of  flights  between  them  allows 
for  the  computation  of  the  Av  required  at  each  bum  point.  Hill’s  equations  (16:80)  were 
used  to  determine  the  total  impulse,  Avt,  required  for  a  particular  maneuver.  This 
parameterization  assumes  the  time,  ti,  at  each  point  is  known;  the  times,  ti^  are 
independent  variables. 

A  complete  circumnavigation  is  defined  by  a  2%  rotation  in  y  from  some  given 
initial  position  (defined  by  a  Yo  on  the  nominal  path)  within  the  required  total  time  of 
flight  (TOP).  The  parameter,  b,  indicates  the  total  number  of  discrete  burn  points  along 
the  circumnavigation  path. 

A  state  vector,  X,  is  composed  of  the  spatial  degrees  of  freedom  for  the  bum 
points’  positions,  and  the  corresponding  times  when  the  deputy  is  located  at  the  bum 
point  positions.  There  are  two  probable  cases  investigated  for  constraining  the  placement 
of  the  burn  points.  First,  the  ‘special  case  ’  requires  all  the  burn  points  to  be  placed  on  the 
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nominal  path.  Therefore,  the  special  case  only  requires  one  degree  of  freedom,  yi,  to 
define  a  bum  point  placement.  Next,  the  ‘general  case  ’  is  defined  by  allowing  the  burn 
points  to  be  placed  anywhere  within  the  constraint  volume.  For  both  cases,  the  burn 
point  timing  is  not  specified,  only  the  sum  of  the  times  as  defined  by  the 
circumnavigation. 

Special  Case  State  Vector. 

The  special  case  state  vector  is  built  from  the  discrete  values  of  yi  and  ti,  defined  by 


X  = 


Yi 

Y,' 


,  for  i  = 


(6) 


where  b  is  the  total  number  of  burn  points.  In  order  to  assure  the  circumnavigation  is 
complete  (i.e.  the  deputy  returns  to  the  initial  point),  the  angle  to  the  final  position,  yb,  is 


computed  as  27t  minus  the  sum  of  the  previous  yi’s.  Similarly,  the  time  of  the  final  point, 
tb,  is  computed  as  TOF  minus  the  sum  of  all  the  previous  t's.  These  values  are  computed 
as 


b-l 

t,=TOF-f^X(j  +  b-l) 

j=l 


(7) 
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General  Case  State  Vector. 


A  general  case,  representing  three  degrees  of  spatial  freedom,  defines  the  burn 
points  placed  within  a  solid  torus  whose  minor  radius  is  defined  by  Pmax-  Any  point  along 
the  actual  flight  path  (and  within  the  solid  torus)  has  a  vector  from  it  to  the  nominal  path 
which  represents  a  deviation  radius,  and  has  a  magnitude  of  p.  This  radius  is  rotated 
about  the  nominal  path  by  an  angle,  e  as  shown  in  Figure  2.  The  angle  e  can  be  rotated 
through  2p  defining  a  circle  about  any  point  on  the  nominal  path.  Rotating  this  circle  by 
?,  creates  a  torus  with  an  inner  radius,  rc  =  ro-p,  and  an  outer  radius,  rt  =  ro  +  p.  The  irmer 
and  outer  radii  will  be  used  in  Chapter  V  below. 


Figure  2.  Sketch  of  Torus  Parameterization 

Figure  2  gives  an  illustration  of  the  torus  parameters  to  define  a  unique  point 
within  the  torus.  A  coordinate  frame  is  fixed  in  the  torus  where  the  p2  direction  is  defined 
by  the  position  of  the  initial  position,  which  can  be  defined  by  the  initial  angle,  yo,  with 
respect  to  the  LVLH  frame.  The  ps  direction  lies  within  the  nominal  path  plane 
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orthogonal  to  p2,  and  pi  completes  the  triad.  The  radius,  p,  is  allowed  to  vary  from  0  to 

Pmax- 

Any  point  in  the  solid  torus  can  be  given  by 
^pPx  +  JpPi  +  ZpPi  =  P  sin(e)Pi  +(ro  +  p  cos(£)cos(Y))p2  +(^o  +  P  cos(£)sin(y))p3  (8) 


where  Xp,  yp,  and  Zp  are  components  in  the  path  coordinate  system  [pi  p2  ps].  Equation  (8) 
is  modified  from  a  general  torus  parameterization.  (8)  The  position  vector  is  expressed  in 
LVLH  coordinates  using  the  nominal  path  rotation  angles  0y  and  0z  discussed  above. 
Thus  the  position  vector,  5r,  is  defined  in  the  LVLH  frame  as 


5F(p,e,Y) 


(p  sin(£))cos0^)cos0^)  +  (r^  +  p  cos^))sin(/)cos0^)sin(0^)-(r^  +pcos$))costS')sin(0^) 
(r^  +pcos$))cos^)cos0J  +  (p  sin(£))cos0^)sin(0J  +  (r^  +  pcos$))sin(/)sin(0^)sin(0J 
(r  +  pcos^))sin^)cos0^)-(psin(£))sin(0^) 


(9) 


where  Jo  is  set  to  zero. 

Now  the  general  case  state  vector  is  defined  using  the  three  degrees  of  freedom: 


Pdi 

P  di-^l 

Si 


A  = 


,  for  /  =  1,-,  b  -1 


(10) 


Y,- 

L 


where  the  final  point  in-plane  angle,?/b,  and  time,  tb,  are  computed  as  in  Eq.  (7),  but  the 


final  point  radius,  p?and  torus  angle,  2,‘hre  allowed  to  vary. 
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Total  Impulse  Computation  (Cost  Function). 


To  compute  each  discrete  individual  Av  requires  knowledge  of  the  individual  time 
of  flight  between  the  and  the  (/+!)'*'  burn  point,  calculated  as 

(11) 

Once  the  position  of  the  current  burn  point,  6r(ti),  along  with  the  time  of  flight, 
Ati,  to  the  next  bum  point,  6r(ti+i),  is  known;  the  required  velocity,  6v'^(  h),  immediately 
following  the  bum  is  computed  as 

=  -[O^CAO  SFCf.O-ara,,,)]  (12) 

The  velocity  just  prior  to  the  bum,  6v'(  ti),  is  calculated  by 

8v-(f,)  =  <I>„(Af,)  .5r(f,_,)  +  <I>„,(Al,)  ■  8v"(?,_,)  (13) 

where  4>vr  and  d>w  are  defined  above  in  Eq.  (4). (16:80)  This  velocity  is  determined  by 
the  location  and  magnitude  of  the  previous  Av.  Eq.  (13)  is  valid  for  all  bum  points  except 
for  the  initial  one.  At  the  initial  burn  point,  5v'(  ti)  is  assumed  to  be  zero  for  all 
calculations. 

Once  the  velocity  just  prior  to  the  bum  point  and  the  velocity  just  after  the 
impulsive  bum  is  known,  the  Av  magnitude  and  direction  is  computed  as 

Av,  =5v  ■"(?,.  )-5v'0,.)  (14) 

Now  that  the  individual  Av  vectors  are  computed,  the  total  required  impulse  can  be 
minimized.  The  total  impulse  is  given  by 
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(15) 


F{X)  =  Av,  =^||Av,. 

/=0 

Singularity  in  Cost  Function. 

The  cost  function  (Eq.  (15))  is  not  continuous  within  the  feasible  region  of  most 
cases.  This  is  a  result  of  using  the  inverse  of  the  Orv  matrix  in  Eq.  (12);  which  is 
computed  as 


-4sin(\|/)) 


<I>rv  = 


—  8  +  8  cos^ )  +  3vj/  sin(\j/ ) 
2n(cos(q/ )  - 1) 

-8  +  8  cos(q/ )  +  3v|/  sin(\|/ ) 


-  2n(cos(q/)  - 1) 

—  8  +  8  cos(^ )  +  sin(\j/ ) 
n  sin(\|/ )  -  1) 

8-8  cos(q/ )  -  3v|/  sin(\|/ ) 


0 


0 


0 

0 

n 

sin(\|/ ) 


(16) 


where  ijr  =  n  t.  This  matrix  contains  a  singularity  when  \|/  is  zero  or  an  integer  multiple  of 


7t.  The  cost  function’s  gradient  becomes  very  steep  in  the  region  near  the  singularity,  and 
the  numerical  optimization  routine  will  not  converge  to  a  solution  if  the  search  routine 
approaches  the  singularity.  Physically,  this  singularity  is  represented  by  burning  between 
two  points  of  a  finite  distance  apart  in  zero  time  which  requires  infinite  Av. 

Several  strategies  are  employed  to  mitigate  the  singularity’s  effects  on  the 
numerical  optimization.  These  strategies  include  initial  guess  inputs  into  the  optimization 
close  to  a  local  minimum,  utilizing  a  low  number  of  burn  points  in  the  initial  guess, 
searching  over  different  potential  numbers  of  burn  points,  and  establishing  bounds  on  the 
time  instances  to  be  £  <  t,.  In  practice,  this  lower  limit  of  the  times,  £,  has  been  set  to 
(1x10'^)  X  TOP  to  prevent  the  routine  from  approaching  too  close  to  the  singularity. 
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These  mitigations  allow  good  performance  of  the  optimization  routine  in  the  vicinity  of  a 
local  minimum. 

Flight  Path  Constraint 

As  mentioned  previously,  the  actual  circumnavigation  flight  path  between  burn 
points  must  be  contained  within  the  torus  and,  therefore,  must  not  deviate  from  the 
nominal  path  by  more  than  pmax-  For  this  reason,  the  deviation  from  the  nominal  path 
must  be  calculated  to  ensure  the  relative  satellite  trajectory  meets  this  constraint.  The 
flight  path  between  burn  points  is  calculated  in  discrete  time  steps  by  propagating  Hill’s 
equations  (16:80)  forward  after  the  burn  is  applied: 


6r  )  =  )  •  6r  (f .  )  +  <!),,  )  •  6v  ^  )  (17) 

where  Afnt  are  intermediate  time  steps  defined  by  dividing  Ati  by  z  time  steps.  A  value  of 
z  =  20  was  chosen  for  all  calculations.  This  value  of  z  provided  adequate  resolution  of 
the  path’s  shape  and  magnitude.  The  intermediate  position  vector,  5r(tint),  is  used  to 
determine  the  flight  path  deviation  represented  as  vector,  pdev,  illustrated  in  Figure  3. 
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Figure  3.  Path  Constraint  Definition  Sketch 
The  ‘path’  frame  [pi  p2  p3]  is  defined  by  the  space- fixed  1-2-3  rotation  by  Jo,  0y, 

and  0z  from  the  LVLH  frame.  Assuming  the  initial  position  always  lies  on  the  nominal 
path,  the  corresponding  transformation  matrix,  R,  from  the  path  frame  to  the  LVLH 
coordinates  is 


R  = 


cos(0j)cos(0j) 

cos(  0^ )  sin(  0J )  sin(Yo )  -cos(  )  sin(  0  ^ ) 
cos(0  ^ )  sin(  0^ )  cos(Yo )  +  sin(Yo )  sin(  0  ^ ) 


cos(0^)sin(0^) 

cos(Yo )  cos(0  J  +  sin(Yo )  sin(  0  ^ )  sin(  0  ^ ) 
-  sin(Yo )  cos(  0^ )  -  cos(Yo )  sin(  0^ )  sin(  0^ ) 


-sin(0^) 

cos(0j)sin(Yo) 

cos(Yo)cos(0^) 


(18) 


This  matrix,  R,  is  then  used  in  the  subsequent  Eqs.  (19),  (20),  and  (21).  Each  intermediate 
point  is  transformed  into  the  path  frame  by 


^Pi 

dp, 


(19) 


The  intermediate  position  vector,  6r(tint),  is  projected  into  the  plane  of  the 
nominal  path  to  calculate  the  position  vector  of  the  nominal  path’s  closest  point.  Since 
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there  is  no  preferred  timing  of  the  intermediate  points  along  the  circumnavigation,  the 


closes  point  is  found  as 


(20) 


The  flight  path  deviation  vector,  Pdev,  is  calculated  by  differencing  the  intermediate 
position  vector  and  the  projected  path  vector  in  the  path  frame: 


Pder  =  -SrCfin,))  (21) 

Finally,  the  flight  path  constraint  is  defined  by  deviation  radius  magnitude,  Pdev,  which 
cannot  exceed  pmax-  The  constraint  is  computed  by  ensuring  the  difference  is  never 
positive: 


P 


dev 


P  max 


<0 


(22) 


Additional  Constraints 

Each  burn  point’s  y  and  t  cannot  exceed  2p  and  TOE  respectively  from  the 
problem  statement  above.  Additionally,  the  values  for  y  and  t  must  be  zero  or  greater 
(positive).  These  values  are  used  to  define  the  upper,  Xu,  and  lower,  X,  limits  on  the 
state  vector. 

A  linear  inequality  constraint  is  needed  to  ensure  Jb  and  tb  are  not  negative.  Erom 
Eq.  (7),  it  can  be  seen  the  sum  of  the  y’s  and  the  sum  of  the  t’s  cannot  be  greater  than  2ti 
or  TOE  respectively.  Physically,  this  would  represent  the  circumnavigation  doubling 
back  on  itself  producing  a  negative  jb,  and  time  flowing  backwards  giving  a  negative  tb. 
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Numerical  Optimization 

A  numerical  optimization  technique  is  employed  to  locate  locally  minimum  fuel 
trajectories.  The  general  form  of  the  optimization  problem  can  be  represented  as 

minF(X) 

XeR" 

(23) 

where  Gfcineq(X)  represents  the  nonlinear  constraints  from  Eq.  (22)  and  linear  inequality 
constraints  from  above.  The  states,  55  and  Xu,  represent  the  lower  and  upper  bound 
constraints  on  the  state  vector.  The  goal  is  to  find  the  state  vector  producing  the 
minimum  value  of  the  cost  function,  F(X)  computed  in  Eq.  (15). 

The  cost  function,  F(X),  is  highly  non-linear.  The  Jmincon’  function  in 
MATLAB’s  Optimization  Toolbox  (6)  was  chosen  to  perform  the  optimization.  This 
routine  is  designed  using  Sequential  Quadratic  Programming  (9:3-26)  which  allows  for 
the  use  of  nonlinear  constraints,  is  appropriate  for  a  single  objective  cost  function,  and 
allows  for  linear  constraints  as  well.  It  is  limited  by  the  fact  the  cost  function  must  be 
continuous  over  the  interval,  and  will  also  attempt  to  minimize  the  maximum  constraint  if 
there  is  no  feasible  solution  (9:3-26). 

MATLAB  Optimization  Routine. 

As  mentioned  above,  MATLAB’s  '/mincon’  function  uses  the  Sequential 
Quadratic  Programming  method  which  is  composed  of  three  main  steps:  updating  the 
Hessian  matrix  of  the  Lagrangian  equation,  solving  the  Quadratic  Programming  sub- 
problem,  and  performing  a  Line- search  and  Merit  function  calculation.  (9: Sec. 3,  26). 
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The  optimization  problem  can  be  reformulated  into  the  Kuhn-Tucker  Equations  (9:  Sec.3, 
26): 

m 

VF(X‘)+£a:.VG„,„(X-))  =  0 

i.i  (24) 

where  Xi*  is  a  Lagrange  multiplier  termed  a  Kuhn-Tucker  point  at  a  unique  state  vector, 
X*,  and  m  is  the  number  of  constraints.  This  equation  essentially  balances  the  gradients 
of  the  active  constraints  and  the  gradient  of  the  cost  function  to  find  a  minimum.  If  the 
minimum  lies  on  the  constraint  boundary,  it  may  not  be  a  true  minimum,  but  the  least 
cost  function  value  along  the  boundary.  The  Kuhn-Tucker  points  can  be  found  by 
solving  the  Lagrangian  equation  (9:  Sec.3,  27): 

m 

L(X,X)  =  F(X)  +  £(X,  •G„,,,(X))  (25) 

(=1 

The  algorithm  starts  with  an  initial  guess  state  vector.  From  the  initial  guess,  a 
Hessian  is  computed  using  finite  difference  calculations.  The  Quadratic  Sub-Problem  is 
then  solved  (9:  Sec.3,  28  and  2:238)  to  determine  the  search  direction.  Once  the  search 
direction  has  been  defined  a  line  search  and  merit  function  are  used  to  determine  the  step 
size  in  order  to  update  the  state  vector.  The  updated  state  vector  is: 

(26) 

where  dk  gives  the  search  direction  and  a  is  the  distance  along  the  search  direction.  Once 
the  step  size  is  determined  the  gradient  of  the  function  is  evaluated  at  the  new  point,  and 
evaluated  against  Eq.  (24).  If  the  convergence  criteria  are  not  met,  the  BEGS  (Broyden- 
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Fletcher- Goldfarb-Shanno)  method  (1:330  and  9:  Sec.3,  30)  is  then  used  to  determine  an 


updated  Hessian  matrix,  and  the  procedure  is  reiterated. 

Practical  MATLAB  Usage. 

MATLAB  uses  finite  differencing  techniques  to  numerically  compute  the  gradient 
and  Hessian  of  the  cost  function  at  any  given  state  vector.  These  finite  differences 
require  bounds  on  the  step  size  for  the  finite  difference.  The  default  was  set  at  1x10'^,  but 
1x10'^  produced  higher  quality  results.  Additionally,  the  tolerance  on  the  state  vector,  the 
cost  function  and  the  constraints  can  be  set  as  well.  Changing  these  tolerances  produced 
significant  differences  in  the  output  of  the  optimization  program.  Setting  all  tolerances 
equal  to  1x10'^  produced  the  most  consistent  results  for  all  the  cases  presented. 

Functionally,  the  state  vector  quantities  were  normalized  for  input  into 
MATLAB ’s  ‘finincon’  routine.  This  gives  a  similar  magnitude  between  the  states  and 
produced  better  convergence.  The  angles  e  and  y  were  normalized  by  In  to  give  values 
between  zero  and  one  ,  the  deviation,  p  was  normalized  by  the  maximum  deviation  radius, 
Pmax,  and  time  was  normalized  by  TOF.  (Angles  of  the  special  case  shown  below  were 
not  normalized.) 

Optimization  Results  Check. 

In  order  to  check  the  optimization  program,  the  cost  function  is  evaluated  in  select 
directions  around  the  area  of  the  minimum  given  by  the  program.  The  goal  of  the  check 
is  to  gain  confidence  in  the  optimization  results  and  understand  what  levels  in  the 
numerical  tolerances  produced  the  best  results.  This  method  is  based  upon  a  subset  of  the 
Weierstrass  Theorem.  (1:83) 
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specifically,  the  cost  function  is  evaluated  at  the  optimized  state  vector  stepped  in 
only  one  state  by  a  range  of  finite  steps  and  does  not  utilize  the  full  set  of  possible  states 
as  required  by  the  Weierstrass  Theorem.  Each  step  creates  a  new  state  vector,  X : 

(27) 


where  m  is  the  step  size,  and  g  is  given  as 


0 

0 


gj{m)  = 


m 

0 


(28) 


where  j  represents  the  state  to  be  stepped.  This  step  vector  is  the  same  length  as  the 
original  state  vector,  X. 

If  the  optimized  state  vector  is  to  produce  a  local  minimum  within  the  defined 
tolerances,  the  new  stepped  value  of  the  cost  function,  F(X),  cannot  be  less  than  the  cost 
function  value  from  the  optimization  output: 

F(X')-F{X)>0  (29) 


This  step  method  described  doesn’t  fully  ensure  a  local  minimum  has  been  found 
because  the  cost  function  could  decrease  in  a  direction  not  orthogonal  to  the  states.  For 
instance,  simultaneously  stepping  states  1  and  2  by  some  amount  produces  another 
unique  state  vector  and  thus  another  unique  cost  function  value.  This  check  describes  a 
necessary  condition  for  a  minimum,  but  doesn’t  describe  a  sufficient  condition. 
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If  given  a  stepped  value  of  the  cost  function  for  which  Eq.  (29)  is  not  satisfied,  the 
state  X  does  not  represent  a  minimum  in  the  cost  function.  However,  the  new  state  vector 
must  be  evaluated  to  determine  if  it  meets  the  constraints  as  well.  If  the  new  state  does 
not  meet  the  constraints,  it  is  not  a  valid  state,  and  therefore  must  be  disregarded. 

Functionally  the  new  stepped  state  is  compared  to  the  constraints.  If  the  new  path 
violates  constraints  an  integer  number  is  added  to  a  variable  called  the  check  sum.  This 
check  sum  is  zero  if  no  constraints  were  violated,  and  greater  than  zero  if  the  constraints 
were  violated.  The  other  constraint’s  violations  (upper  and  lower  limits  of  the  state 
vector  and  the  linear  inequality  constraint  on  the  sum  of  the  times  and  y’s)  were  included 
in  this  check  sum  as  well. 

An  alternate  way  to  verify  a  minimum  is  to  use  the  Kuhn-Tucker  necessary 
conditions  (1:122).  The  Kuhn-Tucker  method  requires  the  calculation  of  the  gradient  of 
the  cost  function  at  the  specific  points.  However,  this  method  was  not  needed  since  the 
step  method  described  above  was  able  to  provide  enough  confidence  in  the  minima  found 
by  the  optimization  routine. 
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III.  Equal  Angle/Equal  Time  Method 


The  behavior  of  the  cost  function  is  initially  evaluated  with  a  very  simple  analytic 
method  in  defining  the  placement  and  timing  of  the  individual  Avi’s.  This  method 
provides  a  baseline  for  comparing  the  optimization’s  or  other  analytic  methods’ 
performance.  The  Equal  Angle/Equal  Time  (EAET)  method  is  defined  by  the  bum  points 
placed  upon  the  nominal  path  in  equal  angular  displacements  along  the  nominal  path  and 
spaced  equal  times  apart.  Each  angular  location  and  time  is  given  by 

271  TOF 

yi=-r’  ^i=—r-  (30) 

b  b 

where  b,  the  number  of  bum  points,  must  be  specified. 

Figure  4  shows  the  behavior  of  the  EAET  with  the  TOF  set  to  0.1  times  the  orbital 
period  of  the  chief  (about  9.25  minutes).  The  chief  is  in  a  circular  orbit  with  an  altitude 
of  400  km  for  all  calculations  and  results  shown  throughout  this  thesis,  unless  otherwise 
stated.  The  nominal  path  is  oriented  in  the  x-y  plane  (0y  =  90°  and  ©z  =  0°)  with  the 
initial  point  rotated  45°  along  the  path  (yo  =  45°).  The  nominal  path’s  radius  is  set  at  50 
m,  with  the  deviation  constraint,  Pmax,  equal  to  10  m.  Two  paths  are  shown  in  Figure  4. 
The  ‘Min  Actual  Path’  is  using  five  burn  points  (b  =  5),  and  the  ‘Infeasible  Path’  is  using 
four  burn  points  (b  =  4).  The  constraint  surface  is  shown  as  the  circular  grid,  representing 
the  toms’s  surface. 
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Figure  4.  EAET  with  Minimum  Eeasible  Number  of  Burn  Points  (b  =  5).  0y  =  9O°,0z  = 

0°  ,yo  =  45°,  TOE  =  0.1,  and  pmax  =  10  m 

Minimum  Number  of  Bum  Points  for  Given  Pmax- 

For  the  EAET  method  there  exists  a  minimum  b  for  the  required  pmax  constraint  to 
be  satisfied.  The  deviation  from  the  nominal  path  is  apparent,  and  depends  upon  the  b 
used  in  Equation  30.  The  number  of  bum  points,  b,  is  stepped  in  integer  increments  from 
2  until  the  minimum  number  of  burn  points  case  meets  the  flight  path  constraint.  In 
Eigure  4,  the  £>  =  4  case  exceeds  the  maximum  deviation. 
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Avt  Versus  Number  of  Burn  Points. 


The  difference  in  magnitude  of  Avt  (also  called  the  change  in  velocity  in  some 
Figures)  between  the  continuous  and  the  discrete  methods  of  circumnavigation  can  be 
quantified  by  Figure  5.  Figure  5  shows  the  variation  of  Avt  as  a  function  of  the  number  of 
burn  points  for  several  path  orientations.  The  circumnavigation  for  each  case  shown  has 
its  nominal  path  rotation  0z  =  0°,  a  radius  of  50  meters  and  a  TOF  =  0.1  times  the  chiefs 
period  (555  s). 


Figure  5.  Avt  Versus  the  Number  of  Bum  Points. 


The  flight  path’s  deviation  from  the  nominal  path  increases  as  the  number  of  bum 


points  decreases,  which  can  be  seen  qualitatively  from  Figure  4.  As  the  number  of  bum 
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points  increases  the  value  of  Avt  approaches  a  value  that  it  would  be  for  a  continuous 
bum  on  the  nominal  path. 

The  flight  path’s  deviation  from  the  nominal  path  is  also  a  function  of  the  number  of 
bum  points.  Therefore,  the  value  of  Pmax  determines  the  minimum  number  of  burn  points 
that  produce  a  trajectory  within  the  toms.  The  minimum  number  of  bum  points  is 
determined  by  trial-and-error  as  above  and  represents  the  minimum  Avt  for  the  EAET 
method. 

The  EAET  method  is  a  very  simple  algorithm  allowing  for  a  quick  determination  of 
the  order  of  the  required  Avt,  but  further  investigations  show  room  to  further  minimize 
Avt.  However,  due  to  its  simplicity,  the  EAET  method  presents  a  good  basis  for 
measuring  the  performance  of  the  optimization  results.  The  Avt  for  each  subsequent 
optimization  result  and  analytical  design  method  is  compared  to  the  EAET  by  computing 
the  percentage  of  savings  from  the  EAET  method  for  the  given  circumnavigation. 

Comparison  with  Continuous  Control  Method. 

The  EAET  method  allows  a  direct  comparison  to  continuous  control  techniques.  (4) 
The  paradigm  used  differs  from  continuous  control  technique  by  the  allowance  of  the 
intermediate  flight  path  between  burn  points  to  vary  off  of  the  defined  nominal  path. 

The  continuous  control  paradigm  ensures  all  points  on  the  circumnavigation  follow  the 
nominal  path  vice  the  discrete  methods  proposed  require  the  path  to  lie  within  a 
constrained  region  about  the  nominal  circular  path.  The  variance  or  deviation  from  the 
nominal  path  allows  for  considerable  savings  in  the  Avt. 
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Varying  the  nominal  path’s  orientation  and  TOF  also  affects  Avt  required.  The 
Avt  varies  similarly  with  TOF,  0y,  and  0z  as  found  in  the  continuous  method  (4:3)  except 
the  magnitude  of  the  Avt  is  less.  For  instance.  Figures  6  and  7  demonstrate  the  variation 
of  Avt  with  rotating  a  50  m  nominal  path  about  0y  while  varying  TOF  from  0.1  to  1.  The 
graphs  show  the  Avt  variation  as  computed  by  the  EAET  method  with  six  bum  points 
(b=6). 

The  shapes  of  the  curves  are  the  same  as  demonstrated  in  the  continuous  control 
method  (4:6)  with  a  few  exceptions.  The  overall  magnitude  of  the  surface  is  less  than  the 
magnitude  presented  in  the  continuous  case,  except  at  the  minimum  points.  Also, 
because  the  initial  velocity  is  set  to  zero  relative  to  the  chief,  a  finite  amount  of  Av  is 
required  to  put  the  circumnavigation  on  a  natural  motion  trajectory  which  occurs  at  0y  = 
30°  and  120°  with  a  TOE  =  1  and  0z  =  0°.  (1 1:7)  For  these  cases,  the  continuous  case  Avt 
would  be  zero  as  can  be  seen  from  Eq.  (2)  and  as  presented  in  Reference  4.  (4:6) 
However,  the  minimums  and  maximums  still  occur  at  the  same  respective  values  of  0y, 
0z,  and  TOF. 
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Figure  6.  Avt  Surface  From  Varying  TOF  and  0y 


!(  lO  -' 


Figure  7.  Avt  Cross-Sections  From  Varying  TOF  and  0y 


This  comparison  is  extrapolated  into  the  results  for  the  optimization  cases,  and 
subsequent  analytical  methods.  For  cases  with  TOF  less  than  0.5  times  the  chiefs  orbital 
period,  the  minimum  Avt  occurs  at  0y  =  90°  and  0z  =  0°.  This  case  will  be  the  primary 
example  for  the  lower  total  impulse  trajectories  calculated  below. 
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rV.  Optimization  Results 


The  EAET  method  has  two  drawbacks:  the  Avt  is  not  optimal  and  the 
intermediate  flight  path  constraints  cannot  be  enforced  except  by  trial- and-error  selection 
of  the  number  of  burn  points.  To  investigate  the  behavior  of  optimal  maneuvers  that 
satisfy  the  path  constraint,  we  use  numerical  optimization.  The  optimization  results  are 
presented  using  the  two  defined  cases:  the  special  case  and  the  general  case. 

Special  Case  Results 

The  EAET  method  is  used  as  the  initial  guess  in  the  optimization  of  the  cost 
function  shown  in  Eq.  (15).  Several  other  types  of  guesses  were  investigated,  but  they 
mainly  involved  random  choices  of  the  states,  and  did  not  provide  lower  cost  function 
values  than  the  EAET  method.  Eigure  8  shows  the  results  from  using  the  EAET  initial 
guess  with  b  =  5,  Pmax  =  10  m  and  a  TOE  =  0.1.  The  nominal  path  is  oriented  in  the  x-y 
plane  to  allow  for  simpler  viewing  of  the  actual  flight  paths;  additionally,  the  x-y  plane 
represents  the  minimum  required  circumnavigation  for  this  TOE  as  shown  in  Eigure  7. 
The  initial  guess’s  (EAET’s)  value  for  the  Avt  is  2.6082  m/s  for  the  path  defined  in  Eigure 
8.  The  final  optimized  Avt  value  is  2.4485  m/s  and  represents  a  4.55%  savings  in  Avt. 
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Figure  8.  Five  Burn  Points,  0y  =  90°,  0z  =  0°,  Yo  =  45°  with  TOF  =  0.1  and  Pmax  =  0.01  km 
One  fundamental  differenee  between  the  guess  and  the  optimized  solution  is  the 
deviation  of  the  intermediate  flight  path  from  the  nominal.  Figure  9  shows  the  magnitude 


the  deviation  from  the  nominal  path,  pdev,  for  the  both  the  initial  guess  and  the  optimized 


solution.  The  zero  points  for  both  lines  in  Figure  9  represent  the  bum  points,  whieh  are 
placed  upon  the  nominal  path.  Note  the  EAET  flight  path  does  not  directly  touch  the 
constraint  surface. 
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Figure  9.  Five  Bum  Points  Deviation,  0y  =  90°,  0z  =  0°,  Jo  =  45°  with  TOF  =  0.1,  and 

Pmax  ~  0.01  km 

The  optimized  intermediate  flight  path  touches  the  constraint  boundary  after  the 
initial  burn  and  the  final  bum  as  seen  in  Figure  9 ;  this  excursion  to  the  constraint 
boundary  skirts  the  inner  radius  of  the  torus  as  seen  in  Figure  8.  The  optimization 
consistently  found  minima  where  the  intermediate  flight  path  skirts  the  constraint 
boundary.  This  skirting  is  characteristic  of  all  the  reasonable  optimal  solutions  computed. 

The  next  step  is  to  investigate  the  effect  of  varying  the  flight  path  constraint,  Pmax. 
Figure  10  shows  the  optimization  for  a  path  with  the  same  orientation  and  TOF  as  Figure 
8,  but  with  Pmax  =  20  m.  The  minimum  number  of  burn  points  for  the  EAET  case  is  now 
four.  The  optimized  Avt  was  determined  as  2.05  m/s  which  represents  a  12.5%  savings 
over  the  EAET  value  of  2.34  m/s. 
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Figure  10.  Four  Bum  Points,  Decrease  Constraint:  0y  =  90°,  0z  =  0°,  Yo  =  45°,  TOF  =  0.1 

and  pmax  =  0.02  km 

The  optimization  routine  placed  the  first  and  last  bum  points  such  that  the  first 
and  last  paths  were  tangential  to  the  inner  constraint  surface.  This  was  a  common  theme 
while  decreasing  the  constraint  (increasing  Pmax)-  Ultimately,  increasing  pmax  decreases 
the  total  number  of  points  possible,  and  if  the  maximum  deviation  is  large  enough  the 
circumnavigation  requires  only  two  bums  which  represents  the  minimum  number  of 
burn  points  for  a  reasonable  circumnavigation. 
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Optimization  Results  Check  Output. 


The  optimization  results  check  as  described  above  is  used  to  determine  whether  or 
not  the  optimization  routine  found  a  local  minimum.  A  variety  of  steps  sizes  were 
utilized  ranging  from  1  x  10'^*^  to  1  x  10  *  step  sizes  in  increments  of  one  order  of 
magnitude.  Each  stepped  state  vector  is  evaluated  to  ensure  the  circumnavigation 
constraints  are  met. 

To  show  whether  the  constraints  were  met  or  exceeded  for  each  new  stepped  state 
vector,  a  Violation  Matrix  is  developed.  The  Matrix  is  calculated  for  each  new  stepped 
state  vector,  X,  and  is  a  row  vector  with  six  columns.  If  the  value  in  the  column  is  zero 
then  the  particular  constraint  corresponding  to  the  column  was  met,  likewise  if  the  value 
is  greater  than  zero  that  constraint  has  been  violated  the  integer  value  number  of  times. 
The  constraint  definitions  are  given  in  Table  1. 


Table  1.  Violation  Matrix  Definition 


Violation  Matrix 
Column  Number 

Constraint 

1 

b-l 

LY;  <  271 

i 

2 

b-l 

i 

3 

Yi  <  27t 

4 

ti<  1 

5 

Yi  and  ti  >  0 

6 
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Table  2  shows  part  of  the  step  check  data  from  the  circumnavigation  and 
optimization  presented  in  Figure  8  above.  The  full  data  are  presented  in  Appendix  A. 

The  top  half  of  Table  2  for  a  step  size  of  1x10'^  produces  changes  in  the  cost  function 
much  lower  than  the  cost  function  tolerance;  in  this  case  the  tolerance  on  the  cost 
function  was  set  at  1x10'^.  The  lower  half  (bolded  numbers)  represent  the  step  size  which 
shows  the  optimization  meets  the  local  minimum  requirements.  If  the  Step  Check 
column  is  negative,  then  the  step  size  and  direction  represent  a  more  minimum  value  for 
the  stepped  cost  function.  However,  since  the  optimization  result  was  constrained  the 
bolded  negative  values  all  exceed  constraints  and  therefore  are  not  feasible  states.  The 
feasible  states  are  all  positive.  This  is  true  for  the  step  size  of  1x10'^,  but  decreasing  the 
step  size  to  1x10’^  represents  a  case  where  the  step  size  produces  a  difference  in  the  new 
cost  function  two  orders  below  the  numerical  limit  of  the  optimization  criteria. 
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Table  2.  Special  Case  Step  Optimization  Check  Example 


Step  Size,  m 

Violation 

Matrix 

Step  Check, 
F(X’)-F(X) 
(km/s) 

Check 

Sum 

1.00E-07 

0000001 

7.47E-12 

1 

1.00E-07 

0000000 

2.93E-11 

0 

1.00E-07 

0000000 

3.02E-11 

0 

1.00E-07 

0000000 

3.07E-11 

0 

1.00E-07 

0000000 

-7.93E-12 

0 

1.00E-07 

0000000 

-1.93E-1 1 

0 

7 

1.00E-07 

0000000 

-1.98E-1 1 

0 

8 

1.00E-07 

0000000 

-1.79E-11 

0 

1 

-1.00E-07 

0000000 

-7.47E-12 

0 

-1.00E-07 

0000000 

-2.93E-11 

0 

-1.00E-07 

0000000 

-3.02E-11 

0 

-1.00E-07 

0000000 

-3.07E-1 1 

0 

-1.00E-07 

0000000 

7.93E-12 

0 

-1.00E-07 

0000000 

1.93E-11 

0 

7 

-1.00E-07 

0000000 

1.98E-11 

0 

8 

-1.00E-07 

0000000 

1.79E-11 

0 

1 

1.00E-06 

0000001 

7.47E-11 

1 

1.00E-06 

0000000 

2.93E-10 

0 

1.00E-06 

0000000 

3.02E-10 

0 

1.00E-06 

0000000 

3.07E-10 

0 

1.00E-06 

0000001 

-7.93E-1 1 

1 

1.00E-06 

0000001 

-1.93E-10 

1 

7 

1.00E-06 

0000001 

-1.98E-10 

1 

8 

1.00E-06 

0000001 

-1.79E-10 

1 

1 

-1.00E-06 

0000001 

-7.47E-1 1 

1 

-1.00E-06 

0000001 

-2.93E-10 

1 

-1.00E-06 

0000001 

-3.02E-10 

1 

-1.00E-06 

0000001 

-3.07E-10 

1 

-1.00E-06 

0000001 

7.94E-11 

1 

-1.00E-06 

0000000 

1.93E-10 

0 

7 

-1.00E-06 

0000000 

1.98E-10 

0 

8 

-1.00E-06 

0000000 

1.79E-10 

0 

Unreasonable  but  Feasible  Solutions. 


Some  of  the  optimization  runs  produced  unreasonable  circumnavigations,  but  still 
met  the  mathematical  constraints.  For  instance,  placing  an  initial  guess  with  only  two 
bum  points  spaced  hy%  radians  for  the  same  circumnavigation  requirements  used  in 


Figure  8,  yields  an  optimized  path  that  doesn’t  circumnavigate  the  chief,  but  still  lies  in 
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the  feasible  space  of  the  cost  function.  The  initial  guess  is  itself  infeasible,  but  the 
resulting  feasible  optimized  path  is  shown  in  Figure  11.  The  optimized  Avt  is  0.84460 
m/s.  Although  the  optimized  path  does  not  actually  circumnavigate  the  chief,  it  does  meet 
the  constraints  as  stated  above.  This  type  of  minimum  was  only  found  when  the  initial 
guess  is  infeasible. 

0.(Sr 

C.CK 

0.02 

I  0 

-D.t3(2  ■ 

-O.M  - 

i<Knn) 

Figure  11.  Two  Burn  Points.  0y=9O,  0z=O,  yo=45  with  TOF  =  0.1  and  pmax=0.01  km 
As  the  initial  conditions  are  modified  by  placing  the  initial  point  ahead  or  behind 

the  chief  along  the  y-axis,  (yo  =  0  or  7t),  a  minimum  can  be  found  by  moving 
infinitesimally  ‘forward’  along  the  path,  and  infinitesimally  ‘back’  to  the  initial  position. 
The  Avt  for  these  ‘circumnavigations’  approach  zero,  which  represents  the  global 
minimum  for  those  conditions.  Of  course,  these  cases  do  not  represent  valid 
circumnavigations.  This  can  be  avoided  by  ensuring  an  adequate  number  of  bum  points 
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are  used  in  the  initial  guess,  that  they  are  sufficiently  spaced  apart,  and  the  initial  guess  is 


feasible. 


Special  Case  Results  Evaluation. 

Several  other  nominal  paths,  TOP  requirements,  and  constraint  boundaries  were 
investigated  using  the  special  case  shown  in  Table  3.  The  EAET  method  was  used  as  the 
initial  guess  for  all  of  the  outputs. 


Table  3.  Special  Case  Representative  Results 


Total  Impulse  (Avt) 

Initial  Conditions 

To  =  50  m, 
a  =  6778  km 

EAET  A  V, 

Mn.  Burn 

Avt 

Snecial  Case  Optimization  Results 

Avt 

(m/s) 

Pts.  (b) 

(m/s) 

%  of  EAET 

13 

0y=6O°,  ©,=30°,  Yo=45°, 
TOF  =  0.1,p„ax=10in 

5 

2.67 

2.55 

s 

o 

C 

o 

4.3% 

Z. 

Cm 

o 

-M 

-M 

a 

’3 

0y=O°,  0,=O°,yo=O°, 

TOF  =  0.1,  p  ma?r=10  m 

5 

3.02 

2.96 

c 

_o 

2.2% 

.S 

O 

’M 

0^9O°,©,=O°,Yo=45°, 

5 

2.61 

2.49 

> 

a 

TOF  =  0.1,  p  max^l^  ^ 

4.5% 

Cm 

o 

0y=6O°,0z=3O°,yo=45°, 

A 

2.39 

2.11 

C 

TOF  =  0.1,  p  max^^^ 

T- 

11.8% 

’-M 

.3 

1 

0y=6O°,0z=3O°,yo=45°, 

2.85 

2.70 

> 

Q. 

TOF  =  0.1,  P  max=8  m 

o 

5.1% 

Cm 

o 

0y=6O°,0z=3O°,yo=45°, 

5 

1.18 

1.13 

c 

TOF  =  0.2,  p  max=10  m 
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The  optimization  results  for  varying  the  orientation  of  the  nominal  path  shows  the 
limit  for  how  well  the  optimization  performs  over  the  EAET  case.  The  best  the 
optimization  performed  was  when  the  nominal  path  is  in  the  x-y  plane,  as  expected  for  a 
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constant  deviation  constraint.  Conversely  the  smallest  difference  between  the  EAET  and 
the  optimization  occurred  when  the  nominal  path  was  not  rotated  at  all,  which 
corresponds  to  the  y-z  plane  circumnavigation.  These  results  are  consistent  with  Figures 
6  and  7  as  well  as  the  continuous  case  results  (4:3-8).  The  variation  show  in  the  y-z  plane 
nominal  path  gives  credibility  to  extrapolating  the  effects  of  nominal  path  orientation 
variation  as  well  as  TOE  variation. 

Comparing  the  total  impulse  value  for  differing  pmax  constraints,  while  keeping 
the  circumnavigation  path  and  TOE  constant,  shows  that  as  pmax  decreases,  the  Avt 
increases.  This  result  is  consistent  with  the  results  implied  from  the  EAET  analytical 
method;  the  minimum  number  of  burn  points  has  to  increase  to  meet  the  flight  path 
constraint  and  as  the  minimum  number  of  burn  points  increases  the  Avt  increases  as  well. 
Another  expected  result  is  the  Avt  scales  directly  with  the  TOE;  shorter  TOE’s  result  in 
greater  A  Vt. 

General  Case  Results 

The  General  Case  allows  the  bum  points  to  vary  off  of  the  nominal  path  with 
three  spatial  degrees  of  freedom,  but  all  other  constraints  apply.  Eigure  12  shows  the 
results  with  0y  =  90°,  0z  =  0°,  jo  =  45°,  TOE  =  0.1,  and  Pmax=0.01  km,  and  the  minimum 
EAET,  b=5,  as  the  initial  guess  (equivalent  to  Figure  8).  The  optimization  routine  found 
a  local  minimum  with  the  intermediate  burn  points  placed  on  the  nominal  path.  The 
optimized  Avt  doubles  the  savings  from  the  optimized  special  case  at  2.3754  m/s 
representing  8.92%  savings.  The  flight  path  touches  the  inner  constraint  radius  at  four 
points  with  the  optimization  routine  minimizing  the  path  length  between  the  third  and 
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fourth  burn  points.  While  all  the  burn  points  still  lie  on  the  nominal  path,  the  endpoint  is 
now  located  on  the  outside  constraint  boundary. 
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Figure  12.  General  Case  Optimization  Results  for  0y=9O,  0z=O,  Yo=45,  TOF  =  0.1 

and  Pmax=0.01  km 


Figure  13.  General  Case  Optimization  Results  Deviation  for  0y=9O,  0z=O,  Yo=45, 

TOF  =  0.1  and  Pmax=0.01  km 
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Figure  12  shows  an  interesting  result;  all  burn  points  still  lie  on  the  nominal  path 
when  they  are  free  to  vary  off  of  it.  This  indicates  the  cost  function  has  a  local  minimum 
when  the  burn  points  are  placed  on  the  nominal  path  and  corresponds  to  the  special  case 
above. 

Varying  Initial  Guess  Radii. 

The  numerical  results  are  sensitive  to  different  initial  guesses.  The  next  approach 
is  to  run  the  optimization  placing  the  initial  guess  bum  points  away  from  the  nominal 
path.  The  guesses  are  defined  by  using  the  EAET  placement  in  y  and  t  along  a  varying 
radius  away  from  the  nominal  path  radius,  but  within  the  constraint  boundary.  Figure  14 
shows  the  placement  of  the  initial  guess  burn  points  on  the  extreme  outer  constraint 
boundary,  at  a  radius  of  60  meters.  The  end  point  of  the  initial  guess  is  set  on  the  outside 
edge  of  the  constraint  boundary.  The  optimization  routine  finds  a  local  minimum  very 
close  to  the  EAET  guess  placed  on  the  nominal  radius  with  the  bum  points  placed  within 
1-3  meters  off  the  nominal  radius. 
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Figure  14.  General  Case  Optimization  with  Guess  on  Outer  Constraint  Boundary. 
0y=9O°,  0z=O°,  Yo=45°  with  TOF  =  0.1  and  pmax=0.01  km 


Figure  15.  Deviation  of  Guess  on  Outer  Constraint  Boundary.  0y=9O°,  0z=O°,  70=45° 

with  TOF  =  0.1  and  pmax=0.01  km 
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Figure  15  shows  the  deviation  of  the  initial  guess  on  the  outside  constraint  boundary  and 
the  optimization.  Again  the  optimization  finds  the  state  vector  such  that  the  intermediate 
flight  paths  skirt  the  inner  constraint  boundary,  while  trying  to  eliminate  a  bum  point  by 
minimizing  the  path  between  the  third  and  fourth  burn  points.  The  Avt  for  this 
optimization  is  2.3775  m/s  or  8.85  %  savings  on  EAET  which  is  approximately  the  same 
as  the  results  for  the  initial  guess  on  the  nominal  path. 

Placing  the  initial  guess  burn  points  on  the  interior  radius  of  the  constraint  surface 
leads  to  non- convergence  in  the  optimization.  The  TOE  of  0.1  causes  every  flight  path 
point  except  for  the  bum  points  to  be  in  the  cost  function’s  infeasible  region.  No  fast 
TOP  was  found  to  converge  on  a  solution  if  all  burn  points  were  placed  on  the  inner 
constraint  radius. 

The  next  step  is  to  progressively  step  the  initial  guess  radius  (constant  radius  upon 
which  the  EAET  points  are  located)  inward  (toward  the  inner  constraint  radius)  from  ro. 
Figure  14  shows  the  results  when  the  5  point  EAET  guess  is  placed  upon  a  radius  of  48.5 
m  (1.5  m  less  than  the  nominal  of  50  m).  The  initial  guess  Avt  is  found  to  be  2.5324  m/s 
and  the  optimized  Avt  is  2.378  m/s.  The  optimized  solution  represents  8.9408%  savings 
on  the  EAET  method  at  the  nominal  radius.  The  key  feature  of  the  new  optimized 
solution  is  the  fact  that  all  intermediate  paths  now  skirt  the  inner  constraint  radius. 

The  Avt  is  approximately  the  same  as  the  previous  two  cases,  indicating  a  region 
in  which  the  objective  function  is  ‘flat’;  resulting  in  a  numerically  sensitive  search  that 
yields  many  different  local  minima.  These  local  minima  are  only  stationary  points  due  to 
numerical  imprecision. 
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Figure  16.  Guess  Radius  of  48.5  m.  0y  =  90°,  0z  =  0°,  jo  =  45°  with  TOF  =  0.1 

and  pmax  =  0.01  km 


Figure  17.  Guess  Radius  of  48.5  m.  Deviation.  0y  =  90°,  0z  =  0°,  Jo  =  45° 
with  TOF  =  0.1  and  pmax  =  0.01  km 

Placing  the  5  burn  point  initial  guess  on  a  radius  tighter  than  the  average  radius  of 
3  m  shown  in  Figure  17  produces  an  infeasible  guess  that  will  not  converge  to  a  solution 
for  these  parameters;  therefore  the  number  of  bum  points  must  be  increased  in  order  to 
find  feasible  initial  guesses  that  will  converge  to  a  solution. 
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Figure  18  represents  the  guess  with  an  inner  radius  one  third  of  the  distance  from 
the  nominal  radius  to  the  inner  constraint  boundary  while  increasing  b  to  six.  The 
optimization  finds  a  minimum  touching  the  inner  constraint  boundary  between  every 
burn  point.  The  resulting  optimized  intermediate  burn  points  are  placed  at  a  nearly 
constant  Pdev  of  6  meters  as  seen  in  the  deviation  of  the  optimization  in  Figure  19. 
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Figure  18.  Initial  Guess  1/3  Less  Nominal  Radius  and  b  =  6.  0y=9O°,  0z=O°,  Yo=45°with 

TOF  =  0.1  and  pmax=0.01  km 
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Figure  19  Deviation  with  Initial  Guess  1/3  Less  Nominal  Radius  and  b  =  6.  0y=9O°, 
0z=O°,  Yo=45°  with  TOF  =  0.1  and  pmax=0.01  km 


The  one  third  radius  guess’s  Avt  is  2.5923  m/s;  whereas  the  optimized  Avt  is  2.3360  m/s. 
This  represents  a  0.61%  savings  and  10.44%  savings  respeetively  eompared  to  the 
minimum  EAET  at  the  nominal  radius.  The  increase  in  savings  from  the  optimized 
solution  shown  in  Eigure  16,  is  counter  to  the  EAET  conclusion  of  the  lower  the  number 
of  burn  points  the  lower  the  Avt. 

The  initial  guess  radius  is  tightened  to  6.5  meters  while  keeping  the  number  of 
burn  points  constant  at  six.  Eigures  20  and  21  illustrate  the  results  from  this  initial  guess. 
Note  the  initial  guess  is  in  the  cost  function’s  infeasible  region,  but  the  optimization 
routine  is  able  to  converge  on  a  minimum.  The  savings  from  the  initial  guess  (compared 
to  EAET)  is  7.23%  at  a  value  of  2.4197  m/s,  whereas  the  savings  from  the  optimization  is 
10.42%  at  a  value  of  2.3365  m/s.  This  shows  there  is  a  unique  radius  for  a  given  number 
of  burn  points  where  all  the  intermediate  flight  paths  are  tangential  to  the  inner  constraint 
radius.  This  fact  will  be  used  to  develop  an  analytical  method  for  determining  the 
placement  of  the  burn  points. 
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Figure  20.  Initial  Guess  Radius  =  6.5  m  and  b  =  6.  0y=9O°,  0z=O°,  Yo=45°  with  TOF  = 

0.1  and  Pmax=0.01  km 


Figure  21.  Deviation  of  Initial  Guess  Radius  =  6.5  m  and  b  =  6.  0y=9O°,  0z=O°,  Yo=45° 

with  TOF  =  0.1  and  Pmax=0.01  km 
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Out  of  Nominal  Plane  Component. 


Most  of  the  burn  points  do  not  have  a  significant  out-of-nominal-plane  placement 
for  most  rotations.  The  maximum  out  of  nominal  plane  placement  of  the  burn  points 
occurs  when  the  nominal  path  is  not  rotated,  but  left  in  the  y-z  plane.  Figure  22  shows  an 
edge  on  view  of  a  circumnavigation  with  rO  =  50  m,  0y  =  0°,  0z  =  0°,  Yo  =  45°,  and  a  TOF 
=  0.1.  The  maximum  deviation  radius  is  set  at  10  m,  but  the  maximum  deviation  out  of 
plane  is  2  m.  The  initial  guess  is  the  EAET  with  a  radius  one  third  less  than  the  nominal 
(the  same  as  presented  in  Figure  22).  The  small  deviation  is  used  as  a  simplifying 
assumption  in  the  development  of  an  analytical  method  below. 
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Figure  22.  y-z  Plane  Circumnavigation 
General  Case  Results  Evaluation. 

The  main  trend  from  the  general  case  optimization  is  the  intermediate  flight  paths 
all  skirt  the  inner  constraint  radius.  These  circumnavigations  represent  the  least  Avt 
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found  for  a  given  number  of  burn  points,  and  the  shortest  path  lengths  between  the  bum 
points. 

There  are  two  additional  trends.  First,  given  a  set  of  rotation  angles,  constraint 
boundary,  and  an  initial  guess  with  a  tighter  radius  the  number  of  bum  points  must 
increase  to  find  a  feasible  solution.  However,  this  increase  in  the  number  of  burn  points 
does  not  increase  Avt  as  indicated  from  the  EAET  method.  Second,  there  appears  to  be 
radius  which  allows  a  balance  between  the  dynamics  of  the  system  and  shortening  the 
path  length  for  a  given  number  of  bum  points.  Additionally,  the  resulting  out-of- 
nominal-plane  component  for  optimized  burn  points  is  small. 

The  initial  guesses  have  so  far  produced  hints  as  to  where  the  minima  are  located, 
but  a  new  method  for  developing  an  initial  guess  is  desired.  A  better  guess  will  allow  for 
a  more  comprehensive  search  for  the  best  circumnavigation  maneuver. 
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V.  Analytical  Design  Method 


Based  on  the  general  case  optimization  results,  an  improved  analytical  design 
method  is  developed  yielding  a  lower  Avt  than  the  EAET  method.  For  fast  TOF  (less 
than  0.3  of  the  chiefs  orbital  period),  it  is  hypothesized  that  a  minimal  total  length  path 
that  skirts  the  inner  constraint  radius  is  optimal  for  a  given  number  of  burn  points.  The 
path  can  be  approximated  as  straight  lines  tangential  to  the  inner  constraint  radius. 
Furthermore,  there  exists  a  unique  radius  from  the  chief  upon  which  the  bum  points  are 
placed.  This  radius  (termed  the  design  radius,  ra)  is  no  more  than  the  nominal  radius  and 
no  less  than  the  nominal  radius  minus  the  maximum  deviation  radius  (the  inner  constraint 
radius).  Trends  in  the  general  case  results  justify  this  hypothesis. 

The  general  case  results  justify  this  hypothesis  based  upon  three  factors.  First,  the 
general  case  optimization  results  indicate  minima  with  all  the  intermediate  flight  points 
touching  the  inner  constraint  radius.  Next,  the  bum  points  lie  within  a  nearly  constant 
radius  greater  than  the  inner  constraint  radius,  but  less  than  the  nominal  radius.  Finally, 
the  optimized  burn  points  lie  within  a  plane  very  near  the  plane  of  the  nominal  path. 

Analytical  Design  Assumptions  and  Simplifications 

In  creation  of  the  design  algorithm,  several  factors  are  assumed  and 
simplifications  made  from  the  general  case  optimization  results.  The  flight  paths  are 
assumed  to  be  straight  lines.  Additionally,  the  algorithm  places  the  ‘designed’  plane  to 
coincide  with  the  nominal  path’s  plane.  The  times  of  flight  between  bum  points  are 
assumed  to  be  the  same  ratio  with  TOF  as  the  bum  point’s  y  angle  from  the  last  bum 
point  to  27U.  The  end  point  for  the  circumnavigation  is  constrained  to  lie  on  the  outside 
constraint  boundary  on  a  line  between  initial  point  and  the  chief.  Finally,  the  y  angles 
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between  the  2"^*  bum  point  and  the  next  to  last  burn  point  are  assumed  to  be  equally 
divided. 

Algorithm  Description  and  Analysis 

All  burn  points  but  the  initial  burn  point  lie  on  a  circle  with  the  design  radius,  ra. 
The  inner  constraint  radius,  rc,  is  computed  by  subtracting  pmax  from  ro.  A  search  range 
for  the  values  of  the  design  radius,  ra,  is  a  set  of  the  radii  from  the  inner  constraint  radius 
to  the  nominal  radius.  A  set  of  burn  points  is  determined  for  each  potential  design  radius 
within  this  range.  Figure  23  shows  the  corresponding  position  and  relevant  angles  for  all 
the  burn  points  for  a  unique  design  radius. 


The  second  burn  point,  bi,  is  determined  by  finding  the  tangent  line  from  the 
initial  point,  bi,  to  the  inner  constraint  radius  in  the  circumnavigation’s  direction,  and 
placing  b2  at  the  farthest  intersection  of  the  tangent  line  and  the  circle  inscribed  by  ra.  In 
order  to  compute  the  y  angle  between  the  two  bum  points,  the  angle  a  must  be  found  by 
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The  next  step  is  to  compute  the  last  bum  point  (in  Figure  23  it  is  be).  The  angle  ttt  is 
computed  exactly  as  Eq.  (31),  but  ro  is  replaced  with  the  radius  to  the  end  point,  rt.  Eq. 
(32)  is  used  again  to  determine  the  final  path  length,  pb,  however  ro  is  replaced  by  ^  in 
Eqs.  (31)  and  (32).  Next,  Eq.  (33)  is  used  to  determine  the  last  angle,  Yb,  replacing  the 
appropriate  variables. 

Once  the  second  and  last  bum  points  have  been  determined,  the  intermediate  bum 
points  are  calculated.  A  tangent  line  to  the  inner  constraint  radius  is  circled  from  the 
second  burn  point  and  each  subsequent  burn  points.  The  next  bum  point  is  placed  at  the 
intersection  of  the  tangent  line  and  the  design  radius.  The  y  value  for  each  bum  point  is 
determined  by  simple  geometry.  This  procedure  is  accomplished  until  the  sum  of 
previous  y  values  exceeds  2ti  -  yb-i.  The  next  to  the  last  angle,  yb-i,  (in  Eigure  23:  ys)  is 
calculated  as 

^b-2 
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The  next  step  is  to  iterate  the  cost  function  through  each  value  of  the  design  radii. 
This  iteration  can  be  visualized  in  Figure  24.  Figure  illustrates  four  design  radii,  special 
attention  should  be  drawn  to  the  two  middle  circles  (specifically  to  the  drawn  area  of 
interest)  and  the  next- to-the- last  path  length.  For  the  greater  of  the  two  circles,  the  next- 
to- the- last  path  is  tangent  to  the  inner  constraint  radius.  However,  decreasing  the  design 
radius  to  the  next  circle  produces  a  next- to-the- last  path  that  is  not  tangent  to  the  inner 
constraint  radius.  There  exist  specific  design  radii  which  produce  a  flight  path  tangent  to 
the  inner  radius  between  each  burn  point  including  the  next  to  the  last  flight  path. 


Figure  24  Sketch  of  Straight  Line  Flight  Paths  with  Multiple  Design  Radii 

Now  that  the  burn  points  are  placed,  the  total  circumnavigation  velocity  change, 
Avt,  is  calculated  for  a  range  of  design  radii.  The  design  radius  producing  the  minimum 
Avt  is  then  chosen.  Figure  25  shows  the  Avt  as  a  function  of  rj  used  to  determine  the  final 
design  radius. 
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Figure  25.  A Vt  vs  rd.  0y  =  90°,  0z  =  0°,  jo  =  45°  with  TOF  =  0.1  and  pmax  =  0.01  km 

The  minimum  for  this  case  occurs  at  a  rd  of  0.04124  km.  The  vertical  dashed 
lines  represent  transitions  from  one  discrete  number  of  burn  points  to  another  in  the 
analytical  design.  The  jumps  in  Avt  in  Figure  25  are  directly  correlated  to  these 
transitions.  The  reason  for  this  behavior  is  seen  in  Figure  24  where  the  pb-i  path  is  only 
tangent  at  unique  values  of  the  design  radius  as  discussed  above. 

An  example  of  the  MATLAB  code  for  this  algorithm  is  located  in  Appendix  B. 

Variation  of  Number  of  Burn  Points  with  Design  Radius. 

The  number  of  burn  points  grows  exponentially  as  the  design  radius  approaches 
the  inner  constraint  radius,  which  can  be  seen  in  Figure  26. 
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Figure  26.  Number  of  Burn  Points  vs.  ra.  0y  =  90°,  0z  =  0°,  Jo  =  45°  with  TOF  =  0.1  and 

Pmax  ~  0.01  km 

The  minimum  Avt  does  not  eorrelate  to  the  minimum  number  of  burn  points.  This 
result  seemingly  eontradiets  the  conclusion  above  that  the  lower  the  number  of  burn 
points  the  lower  Avt.  The  number  of  bum  points  is  increased  in  order  to  shorten  the  flight 
path  to  skirt  the  inner  constraint  radius.  Eventually,  though,  increasing  the  number  of 
bum  points  (tightening  the  design  radius)  no  longer  decreases  the  Avt  but  increases  it. 

End  Point  Position  Selection. 

The  decision  to  nominally  place  the  end  point  on  the  outer  constraint  boundary 
was  determined  empirically  by  investigating  many  general  case  solutions,  and 
experimentally  by  varying  the  end  point  for  a  specific  circumnavigation  maneuvers 
analytical  design.  The  analytical  design  algorithm  is  modified  to  change  the  position  of 
the  end  point  for  the  maneuver.  This  is  accomplished  by  varying  the  value  of  p  from  the 
inner  toms  radius  to  the  outer  toms  radius  and  computing  the  total  impulse  for  each 
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design.  Figure  27  shows  the  impact  on  the  total  impulse  with  varying  the  radius  to  the 
end  point  for  the  circumnavigation  maneuver  defined  by  the  parameters:  0y  =  90°,  0z  =  0°, 
Yo  =  45°  with  TOF  =  0.1  and  pmax  =  0.01  km.  The  maximum  value  for  Avt  occurs  on  the 
inner  constraint  radius,  and  the  minimum  value  occurs  at  the  outer  constraint  radius. 


-3 

X  10 


Figure  27.  End  Point  Variation  in  the  Analytical  Design  Method,  Pmax  =  0.01  km. 

There  are  circumnavigations  for  which  placing  the  endpoint  on  the  outer 
constraint  boundary  does  not  represent  the  minimum  analytical  design  choice.  Figure  28 
shows  the  results  of  varying  the  endpoint  for  the  same  circumnavigation,  but  with  an 
increased  maximum  deviation  radius:  pmax  =  0.02  km. 
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Figure  28.  End  Point  Variation  in  the  Analytical  Design  Method,  Pmax  =  0.02  km. 

This  shows  the  optimal  placement  of  the  endpoint  near  the  nominal  radius; 
however  the  difference  between  placing  the  endpoint  on  its  optimal  position  and  the  outer 
constraint  boundary  is  not  great.  Additionally,  the  large  increase  in  the  maximum 
deviation  radius  for  this  circumnavigation  may  not  be  realistic  for  operational  constraints. 
To  simplify  the  algorithm,  the  end  point  is  assumed  to  lie  on  the  outer  constraint 
boundary  for  the  rest  of  the  results  presented.  The  optimal  end  point  placement  can  be 
used  if  the  maximum  deviation  radius  becomes  a  significant  fraction  of  the  nominal 
radius,  however  for  all  examples  in  this  thesis  the  end  point  has  been  placed  on  the  outer 
constraint  boundary. 
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Analytical  Design  Method  Performance 


The  designed  algorithm  for  the  initial  conditions  is  input  as  the  initial  guess  into 
the  optimization  routine.  This  allows  for  significant  savings  compared  to  the  EAET  case, 
but  also  validates  the  use  of  the  algorithm  as  a  good  approximation  near  a  local 
minimum. 

Optimization  Results  with  Design  as  Initial  Guess. 

Figure  29  shows  the  optimization  results  using  the  design  algorithm  guess  for 
ro  =  50  m,  0y  =  60°,  0z  =  30°,  yo  =  45°  with  TOF  =  0.1  and  Pmax=0.01  km.  Figure  30 
shows  the  deviation  of  both  the  design  and  the  optimized  solution. 


•  Bvri  Pcints 

—  Pith 

-I-  OptimiiMl  Bum  Poinfe 

—  Optinmzfrd  Palh 

—  Nominal  Path 


Figures  29.  Design  and  Optimized  Circumnavigation  Paths,  ro  =  50  m,  0y  =  60°, 
0z  =  30°,  Yo  =  45°  with  TOF  =  0.1  and  pmax  =  0.01  km 
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Figures  30.  Design  and  Optimized  Circumnavigation  Path  Deviations,  ro  =  50  m,  0y  = 
60°,  0z  =  30°,  Yo  =  45°  with  TOF  =  0.1  and  Pmax  =  0.01  km 

The  two  paths  are  very  close  to  each  other,  with  the  optimized  path  touching  the 
inner  constraint  radius  at  most  points.  The  design  actual  path  does  not  actually  touch  the 
constraint  radius  because  the  actual  flight  paths  between  bum  points  are  curved  as 
opposed  to  straight  lines.  The  Avt  for  the  design  path  is  2.3928  m/s  with  the  optimized 
Avt  equal  to  2.2942  m/s  which  represent  savings  of  10.30%  and  14.00%  respectively. 
Savings  realized  by  using  the  optimized  state  vector  as  opposed  to  the  design  state  vector 
are  only  0.0985  m/s,  or  3%  on  EAET.  The  computational  requirements  for  the 
optimization  may  not  justify  this  small  increase. 

Variation  of  the  Inner  Constraint  Radius  (Maximum  Deviation  Radius). 

As  seen  in  the  Special  Case  results  varying  the  maximum  deviation  radius  has  an 
impact  on  the  placement  of  the  bum  points  for  the  optimized  solution.  For  the  analytic 
method  the  number  of  bum  points  increases  as  expected  (subsequently  the  Avt  increases) 
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as  the  maximum  deviation  radius  decreases.  The  inverse  relationship  is  true  as  well. 
Figures  3 1  and  32  show  the  circumnavigation  paths  and  deviations  respectively  for  the 
same  circumnavigation  shown  in  Figure  29,  but  with  pmax  increased  to  20  m.  The  Avt  for 
the  design  path  is  1.84  m/s  with  the  optimized  Avt  equal  to  1.72  m/s  which  represent 
savings  of  21.25%  and  28.33%  respectively  over  an  EAET  method  using  4  bum  points 
and  a  Avt  =  2.40  m/s. 


O  Design  Bum  Points 

-  Design  Path 

+  Optimized  Bum  Points 
—  Optimized  Path 
- Nominal  Path 


Figure  31.  Design  and  Optimized  Circumnavigation  Paths — Larger  pmax- 
ro  =  50  m,  0y  =  60°,  0z  =  30°,  yo  =  45°  with  TOF  =  0.1  and  pmax  =  0.02  km 
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Figure  32.  Design  and  Optimized  Circumnavigation  Deviations — Larger  Pmax- 
ro  =  50  m,  0y  =  60°,  ©z  =  30°,  Jo  =  45°  with  TOF  =  0.1  and  Pmax  =  0.02  km 

The  design  radius  for  this  case  is  farther  from  the  inner  constraint  radius  at  32.9  m 
(2.9  m  away  from  the  inner  constraint  radius).  The  major  increase  in  savings  is  a  result 
of  the  total  path  length  for  the  circumnavigation  being  considerably  shorter  (tighter  radius) 
than  the  EAET  method,  even  with  reducing  the  number  of  burn  points.  The  amount  of 
decrease  in  Avt  for  the  EAET  due  to  the  decrease  in  the  number  of  bum  points  (compared 
to  the  previous  one)  does  not  compensate  for  the  better  results  found  by  the  design  and 
optimization.  Overall,  comparing  the  actual  values  of  Avt,  the  general  observation  is 
made:  increasing  the  pmax  results  in  a  decrease  in  Avt,  consistent  with  all  other  cases. 
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Variation  of  the  Nominal  Radius. 


Increasing  the  nominal  radius  shows  the  performance  of  the  algorithm  for  another 
circumnavigation,  and  the  percentage  of  savings  over  the  EAET  method  stays  relatively 
constant.  Figures  33  and  34  show  the  results  of  increasing  the  nominal  radius  to  100m 
while  making  keeping  the  percentage  (of  the  nominal  radius)  of  the  maximum  deviation 
radius  constant  at  20%  (20  m  for  100  m).  The  orientation  of  the  nominal  path  and  TOE 
are  the  same  as  the  circumnavigation  shown  in  Figures  29  and  30  (ro  =  100  m,  0y  =  60°, 
0z  =  30°,  Yo  =  45°  with  TOE  =  0.1  and  pmax  =  0.02  km).  The  Avt  for  the  design  path  is 
4.78  m/s  with  the  optimized  Avt  equal  to  4.59  m/s  which  represent  savings  of  10.30%  and 
13.90%  respectively. 


Figures  33.  Design  and  Optimized  Circumnavigation  Paths — Larger  ro. 
ro  =  100  m,  0y  =  60°,  0z  =  30°,  Yo  =  45°  with  TOF  =  0.1  and  pmax  =  0.02  km 
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Figures  34.  Design  and  Optimized  Circumnavigation  Path  Deviations —  Larger  ro. 
ro  =  100  m,  0y  =  60°,  0z  =  30°,  yo  =  45°  with  TOF  =  0.1  and  pmax  =  0.02  km 

There  are  two  effects  of  note  from  this  example.  First,  the  percentage  of  savings 

from  the  EAET  stays  relatively  constant.  Keeping  the  same  ratio  of  inner  constraint 

radius  to  the  nominal  radius  allows  for  relatively  constant  savings,  which  implies  for  any 

given  orientation  the  results  can  be  scaled  at  a  particular  TOF. 

Second,  the  number  of  bum  points  stays  relatively  constant  for  both  the  EAET  and 
the  design/optimization.  As  the  nominal  radius  is  increased  for  a  constant  TOF,  the  time 
to  travel  between  subsequent  bum  points  decreases  with  increasing  distance.  Essentially, 
this  indicates  the  paths  between  the  burn  points  better  approximate  a  line,  and  the 
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ininimum  number  of  points  allowed  with  a  given  inner  constraint  radius  stays  nearly 
constant. 


Circumnavigation  Requirements  Comparison. 

The  results  from  a  representative  set  of  initial  conditions  are  shown  in  Table  4. 
The  analytic  design  outperforms  the  minimum  burn  EAET  case,  and  the  optimization 
results  using  the  analytical  design  method  as  the  initial  guess  always  produces  the 
minimum  circumnavigation  Avt  found  from  any  guess.  Another  benefit  to  the  analytical 
design  as  an  initial  guess  is  a  very  high  probability  of  convergence  to  a  minimum. 


Table  4.  Analytic  Design  Algorithm  Performance 


Initial  Conditions 
ro  =  50  m, 
a  =6778  km 

Total  Impulse  (Avt) 

EAETAvr 

Analytical  Design  Algorithm 

Min. 

Burn 

Pts.  (b) 

^  .  Avt  Optimization 

Design 

Burn  rd  / 

Points  rmt  -  - 

X  OlIllS  vlU/ 

(b)  %  of  EAET 

Variation  of  Nominal 
Path  Orientation 

0y=6O°,  ©z=30°,  Yo=45°, 
TOF  =  0.1,  p„ax=10  m 

5  2.67 

2.39  2.29 

11  11  ^ 

10.3%  14.0% 

0y=O°,  0,=O°,Yo=O°, 

TOF  =  0.1,  p  niax=10  m 

5  3.02 

2.79  2.71 

11  /1 1  ^ 

7.8%  10.5% 

©y=9O°,0,=O°,Yo=45°, 
TOF  =  0.1,  Pmax=10  m 

5  2.61 

2.37  2.24 

11  41  ^ 

10.8%  14.3% 

Variation  of 

Pmax 

0y=6O°,  ©,=30°,  Yo=45°, 
TOF  =  0.1,  p„ax=20  m 

4  2.39 

1.84  1.73 

13  3^^  9 

23.2%  27.8% 

0y=6O°,0z=3O°,yo=45°, 
TOF  =  0.1,  p  max=8  m 

6  2.84 

2.54  2.44 

13  1^  9 

10.9%  14.1% 

Variation  of 
TOF 

0y=6O°,0z=3O°,yo=45°, 
TOF  =  0.2,  p  max^lO 

5  1.18 

1.05  0.99 

10  ^ 

11.3%  16.6% 

0y=6O°,0z=3O°,yo=45°, 
TOF  =  0.05,  pmax=10  m 

5  5.67 

5.14  4.97 

1  ^  /lO  ^  - 

iJ  T-VJ.O 

9.4%  12.3% 
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Table  4  also  presents  other  indications  of  the  cost  function’s  nature,  specifically; 
there  is  a  strong  correlation  between  the  savings  and  the  inner  constraint  radius  as  well  as 
the  TOP.  As  Figure  35  shows,  the  percent  savings  between  the  EAET  method  and  the 
Analytical  Design  Method  decreases  with  increasing  TOE.  As  the  TOE  increases  the 
EAET  becomes  a  better  approximation.  The  local  minimums  produced  using  the  design 
algorithm  as  the  initial  guess  are  the  lowest,  valid  circumnavigations  found  of  all  the 
initial  conditions  investigated. 

The  discontinuous  steps  in  the  EAET  curve  represent  changes  in  the  minimum 
number  of  burn  points.  The  EAET  does  not  necessarily  minimize  the  path  length  by 
having  the  intermediate  flight  path  tangent  to  the  inner  constraint  radius,  which  produces 
the  major  difference  between  the  two  curves.  The  analytical  method  has  the  advantage  of 
producing  a  relatively  smooth  line  in  Figure  35  which  eliminates  inefficiencies  created  by 
the  discretized  nature  of  the  EAET  method. 


Eigure  35.  Avt  versus  rc  with  varying  TOF. 
ro  =  50  m,  0y=6O°,  0z=3O°,  Yo=45°,  pmax  =  0.01  km. 
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Design  Radius  Versus  Inner  Constraint  Radius  (Maximum  Deviation  Radius). 


The  analytical  design  method  determines  a  design  radius  based  upon  the 
dynamics  and  the  inner  constraint  radius.  It  is  informative  to  investigate  how  the  design 
radius  varies  as  a  function  of  the  inner  constraint  radius  for  a  particular  circumnavigation 
which  is  presented  in  Figure  36. 


Figure  36.  Design  Radius  Versus  Inner  Constraint  Radius,  ro  =  50  m,  0y=6O°,  0z=3O°, 

Yo=45°,  TOF  =  0.1 


The  dashed  straight  line  is  where  the  design  radius  would  be  equal  to  the  inner  constraint 
radius,  and  the  starred  line  represeits  the  computed  design  radius  for  each  inner 
constraint  radius.  The  design  radii  approach  the  inner  constraint  radii  as  the  inner 
constraint  radius  shrinks  with  the  greatest  difference  at  the  30  m  point.  The  greatest 
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difference  is  approximately  3  m.  There  does  not  appear  to  be  a  direct  correlation  between 
the  differences  in  the  design  radii  and  the  inner  constraint  radii  and  the  total  impulse 
shown  in  Figure  35.  The  roughness  of  the  design  radius  line  is  a  function  of  the 
increment  size  of  the  individual  design  radii  used  in  the  algorithm;  the  smaller  the 
increment  size  the  smoother  the  line. 

Initial  Conditions  Considerations. 

The  analytical  design  method  can  be  easily  generalized  to  the  problem  where  the 
initial  position  does  not  necessarily  lie  on  the  ‘nominal’  path.  This  can  be  accomplished 
by  replacingthe  nominal  radius  with  a  ‘pseudo-nominal’  radius  in  the  algorithm,  and 
placing  the  inner  constraint  radius  at  the  minimum  allowable  distance  to  the  chief.  For 
this  problem,  the  end  point  placement  would  need  to  be  taken  into  account.  If  the 
distance  from  the  initial  point  to  the  inner  constraint  radius  (i.e.  maximum  deviation 
radius)  is  large  (greater  than  0.2  ro),  then  the  position  of  the  end  point  should  be  added  to 
the  algorithm  as  shown  above. 

In  general,  the  initial  velocity  of  the  deputy  with  respect  to  the  chief  will  not  be 
necessarily  zero.  The  plane  of  the  constraint  torus  may  be  chosen  (if  a  free  parameter)  to 
better  the  deputy’s  initial  relative  velocity  vector.  The  magnitude  and  direction  of  the 
initial  velocity  can  have  a  significant  affect  on  the  Avt  by  directly  adding  or  subtracting  to 
the  results  presented. 
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VI.  Conclusions 


Two  design  approaches  have  been  developed  to  allow  a  satellite  to 
circumnavigate  a  chief  satellite  or  Resident  Space  Object  via  impulsive  thrusting.  The 
first  method  located  burn  points  about  a  nominal  path  in  equal  angles  and  equal  times. 
When  considering  the  flight  path  constraints,  there  exists  a  specific  minimum  number  of 
burn  points  for  which  the  equal-angle/equal-time  flight  path  will  be  feasible.  This 
minimum  number  of  bum  points  yields  the  minimum  total  fuel  expenditure  for  this 
method.  This  method  is  useful  due  to  its  inherent  simplicity,  and  scales  to  the  continuous 
case  as  the  number  of  burn  points  goes  to  infinity. 

The  second  circumnavigation  analytic  approach  was  to  choose  a  design  radius  and 
locate  the  burns  along  a  circle  of  that  radius,  with  the  trajectory  between  the  burns 
intersecting  a  circle  of  minimum  approach  distance  at  a  tangent  point.  This  method 
always  yielded  lower  fuel  expenditure  than  the  equal-angle/equal-time  (EAET)  method, 
and  increased  the  savings  over  the  EAET  method  as  the  total  circumnavigation  time  of 
flight  decreased.  The  analytical  design  method  produced  a  relatively  simple  algorithm 
that  could  be  used  for  on-board  autonomous  operations. 

The  two  described  methods  were  employed  as  initial  guesses  for  a  numerical 
optimization  routine  to  find  a  minimum- fuel  solution.  In  most  cases,  the  local  minimum 
was  fairly  close  to  the  analytical  design  method  in  both  trajectory  and  fuel  expenditure. 
Thus,  the  analytical  design  approach  generally  always  provides  the  best  initial  guess 
toward  a  minimum- fuel  solution. 


67 


Appendix  A.  Special  Case  Output  Data  with  Optimization  Check 


Special  Case  Data  Output : 


Optimized  Total  Delta  Vee: 
Initial  Total  Delta  Vee: 
Inputs: 

thetayd  =  90  deg 

thetayzd=  0  deg 

gamma0d=  45  deg 


b=  5 

ro_m=  50  km 

a=  6778  km 

Tot_TOF=  0.100 

r_dev  =  0.0100  km 

TolFun  =  l.Oe-009 

TolCon=  l.Oe-009 

TolX=  l.Oe-009 

DiffMaxChange  =  l.Oe-001 
DiffMinChange  =  l.Oe-009 


2.48949 127357e-003 
2.608 17455 142e-003 


Eigenvalues  of  Hessian: 


Gradient: 

0.00007461 

0.00029368 

0.00030405 

0.00035531 

-0.00006748 

-0.00017411 

-0.00018238 

-0.00007135 


0.00038337 

0.01170769 

0.06065438 

0.12869209 

0.74177054 

0.91848107 

1.03272863 

2.00804976 


Norm  of  Gradient:  0.00061946 


Initial  Guess  State  Vector: 

1.25664 

1.25664 

1.25664 

1.25664 
0.20000 
0.20000 
0.20000 
0.20000 


Optimized  State  Vector: 

1.43005 

1.15891 

1.14946 

1.14303 

0.25428 

0.19841 

0.18950 

0.18143 
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X 

g(u) 

Violation  Mat  F(X') 

F(X’)-F(X) 

Ck  Sum 

1 

+1.0e-009 

0000000 

2.48949 13e-003 

+7.47245 15e-014 

0 

2 

+1.0e-009 

0000000 

2.48949 13e-003 

+2.9334521e-013 

0 

3 

+1.0e-009 

0000000 

2.48949 13e-003 

+3.0243429e-013 

0 

4 

+1.0e-009 

0000000 

2.48949 13e-003 

+3.0654038e-013 

0 

5 

+1.0e-009 

0000000 

2.48949 13e-003 

-7.9358395e-014 

0 

6 

+1.0e-009 

0000000 

2.48949 13e-003 

-1.9296023e-013 

0 

7 

+1.0e-009 

0000000 

2.48949 13e-003 

-1.9792024e-013 

0 

8 

+1.0e-009 

0000000 

2.48949 13e-003 

-1.7930579e-013 

0 

1 

-l.Oe-009 

0000000 

2.48949 13e-003 

-7.47245 15e-014 

0 

2 

-l.Oe-009 

0000000 

2.48949 13e-003 

-2.9334521e-013 

0 

3 

-l.Oe-009 

0000000 

2.48949 13e-003 

-3.0243473e-013 

0 

4 

-l.Oe-009 

0000000 

2.48949 13e-003 

-3.0654082e-013 

0 

5 

-l.Oe-009 

0000000 

2.48949 13e-003 

+7.9305052e-014 

0 

6 

-l.Oe-009 

0000000 

2.48949 13e-003 

+1.9295633e-013 

0 

7 

-l.Oe-009 

0000000 

2.48949 13e-003 

+1.9785952e-013 

0 

8 

-l.Oe-009 

0000000 

2.48949 13e-003 

+1.7924507e-013 

0 

H: 

1 

+1.0e-008 

0000000 

2.48949 13e-003 

+7.4724732e-013 

0 

2 

+1.0e-008 

0000000 

2.48949 13e-003 

+2.9334547e-012 

0 

3 

+1.0e-008 

0000000 

2.48949 13e-003 

+3.0243473e-012 

0 

4 

+1.0e-008 

0000000 

2.48949 13e-003 

+3.065405 le-012 

0 

5 

+1.0e-008 

0000000 

2.48949 13e-003 

-7.9335540e-013 

0 

6 

+1.0e-008 

0000000 

2.48949 13e-003 

-1.9296049e-012 

0 

7 

+1.0e-008 

0000000 

2.48949 13e-003 

-1.9789764e-012 

0 

8 

+1.0e-008 

0000000 

2.48949 13e-003 

-1.7926841e-012 

0 

1 

-l.Oe-008 

0000000 

2.48949 13e-003 

-7.47244716-013 

0 

2 

-l.Oe-008 

0000000 

2.48949 13e-003 

-2.93345256-012 

0 

3 

-l.Oe-008 

0000000 

2.48949 13e-003 

-3.02434556-012 

0 

4 

-l.Oe-008 

0000000 

2.48949 13e-003 

-3.06540386-012 

0 

5 

-l.Oe-008 

0000000 

2.48949 13e-003 

+7.93350196-013 

0 

6 

-l.Oe-008 

0000000 

2.48949 13e-003 

+1.92960196-012 

0 

7 

-l.Oe-008 

0000000 

2.48949 13e-003 

+1.97894656-012 

0 

8 

-l.Oe-008 

0000000 

2.48949 13e-003 

+1.79265636-012 

0 
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X 

g(u) 

Violation  Mat  F(X') 

F(X’)-F(X) 

Ck  Sum 

1 

+1.0e-007 

0000001 

2.48949 13e-003 

-1-7.47245976-012 

1 

2 

+1.0e-007 

0000000 

2.48949 13e-003 

-h2.9334540e-011 

0 

3 

+1.0e-007 

0000000 

2.48949 13e-003 

-h3.0243468e-011 

0 

4 

+1.0e-007 

0000000 

2.48949 13e-003 

-h3.0654051e-011 

0 

5 

+1.0e-007 

0000000 

2.48949 13e-003 

-7.9332938e-012 

0 

6 

+1.0e-007 

0000000 

2.48949 13e-003 

-1.9295616e-011 

0 

7 

+1.0e-007 

0000000 

2.48949 13e-003 

-1.9789163e-011 

0 

8 

+1.0e-007 

0000000 

2.48949 13e-003 

-1.7926180e-011 

0 

1 

-l.Oe-007 

0000000 

2.48949 13e-003 

-7.4724571e-012 

0 

2 

-l.Oe-007 

0000000 

2.48949 12e-003 

-2.9334531e-011 

0 

3 

-l.Oe-007 

0000000 

2.48949 12e-003 

-3.0243460e-011 

0 

4 

-l.Oe-007 

0000000 

2.48949 12e-003 

-3.0654041e-011 

0 

5 

-l.Oe-007 

0000000 

2.48949 13e-003 

-i-7.9339729e-012 

0 

6 

-l.Oe-007 

0000000 

2.48949 13e-003 

-hl.9296525e-011 

0 

7 

-l.Oe-007 

0000000 

2.48949 13e-003 

-hl.9790109e-011 

0 

8 

-l.Oe-007 

0000000 

2.48949 13e-003 

-1-1.79273536-011 

0 

>1:  >1:  >1:  >1:  >1:  >1:  >1: 

1 

+1.0e-006 

0000001 

2.48949 13e-003 

-1-7.47246616-011 

1 

2 

+1.0e-006 

0000000 

2.48949 16e-003 

-1-2.93345726-010 

0 

3 

+1.0e-006 

0000000 

2.48949 16e-003 

-1-3.02435046-010 

0 

4 

+1.0e-006 

0000000 

2.48949 16e-003 

-h3.0654096e-010 

0 

5 

+1.0e-006 

0000001 

2.48949 12e-003 

-7.93019806-011 

1 

6 

+1.0e-006 

0000001 

2.489491  le-003 

-1.92916316-010 

1 

7 

+1.0e-006 

0000001 

2.489491  le-003 

-1.97848546-010 

1 

8 

+1.0e-006 

0000001 

2.489491  le-003 

-1.79207186-010 

1 

1 

-l.Oe-006 

0000001 

2.48949 12e-003 

-7.47245036-011 

1 

2 

-l.Oe-006 

0000001 

2.4894910e-003 

-2.93344996-010 

1 

3 

-l.Oe-006 

0000001 

2.4894910e-003 

-3.02434246-010 

1 

4 

-l.Oe-006 

0000001 

2.4894910e-003 

-3.06539966-010 

1 

5 

-l.Oe-006 

0000001 

2.4894914e-003 

-1-7.93708276-011 

1 

6 

-l.Oe-006 

0000000 

2.48949 15e-003 

-hl.9300532e-010 

0 

7 

-l.Oe-006 

0000000 

2.48949 15e-003 

-1-1.97944196-010 

0 

8 

-l.Oe-006 

0000000 

2.48949 15e-003 

-1-1.79328336-010 

0 
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X 

g(u) 

Violation  Mat  F(X') 

F(X’)-F(X) 

Ck  Sum 

1 

+1.0e-005 

0000001 

2.4894920e-003 

-^7.47253636-010 

1 

2 

+1.0e-005 

0000000 

2.48949426-003 

-^2.93349026-009 

0 

3 

+1.0e-005 

0000000 

2.48949436-003 

-^3.02438626-009 

0 

4 

+1.0e-005 

0000000 

2.48949436-003 

-h3. 06545476-009 

0 

5 

+1.0e-005 

0000001 

2.48949056-003 

-7.89921356-010 

1 

6 

+1.0e-005 

0000001 

2.48948936-003 

-1.92515966-009 

1 

7 

+1.0e-005 

0000001 

2.48948936-003 

-1.97417916-009 

1 

8 

+1.0e-005 

0000001 

2.48948956-003 

-1.78661836-009 

1 

1 

-l.Oe-005 

0000001 

2.48949056-003 

-7.47238026-010 

1 

2 

-l.Oe-005 

0000001 

2.48948836-003 

-2.93341696-009 

1 

3 

-l.Oe-005 

0000001 

2.48948826-003 

-3.02430666-009 

1 

4 

-l.Oe-005 

0000001 

2.48948826-003 

-3.06535446-009 

1 

5 

-l.Oe-005 

0000001 

2.48949216-003 

-^7.96806646-010 

1 

6 

-l.Oe-005 

0000000 

2.48949326-003 

-hi. 93405656-009 

0 

7 

-l.Oe-005 

0000000 

2.48949336-003 

-hi. 98374836-009 

0 

8 

-l.Oe-005 

0000000 

2.48949316-003 

-hi. 79873686-009 

0 

>1:  >1:  >1:  >1:  >1:  >1:  >1: 

1 

+1.0e-004 

0000001 

2.48949876-003 

-h7.47323886-009 

1 

2 

+1.0e-004 

0000000 

2.48952066-003 

-h2.93382046-008 

0 

3 

+1.0e-004 

0000000 

2.48952156-003 

-h3. 02474456-008 

0 

4 

+1.0e-004 

0000000 

2.48952196-003 

-h3. 06590606-008 

0 

5 

+1.0e-004 

0000001 

2.48948376-003 

-7.58930826-009 

1 

6 

+1.0e-004 

0000001 

2.48947246-003 

-1.88512486-008 

1 

7 

+1.0e-004 

0000001 

2.48947206-003 

-1.93112136-008 

1 

8 

+1.0e-004 

0000001 

2.48947406-003 

-1.73209556-008 

1 

1 

-l.Oe-004 

0000001 

2.48948386-003 

-7.47167776-009 

1 

2 

-l.Oe-004 

0000001 

2.48946196-003 

-2.93308666-008 

1 

3 

-l.Oe-004 

0000001 

2.48946106-003 

-3.02394826-008 

1 

4 

-l.Oe-004 

0000001 

2.48946066-003 

-3.06490316-008 

1 

5 

-l.Oe-004 

0000001 

2.48949966-003 

-h8.27784066-009 

1 

6 

-l.Oe-004 

0000000 

2.48951106-003 

-hi. 97409386-008 

0 

7 

-l.Oe-004 

0000000 

2.48951156-003 

-h2.0268 1436-008 

0 

8 

-l.Oe-004 

0000000 

2.48950986-003 

-hi. 85328056-008 

0 
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Appendix  B.  MATLAB  *  Code  for  Analytical  Design  Method^ 


Main  Program:  Analytical  Design  and  Optimization 

%  Analytical  Design  Algorithm  for  Fast  Satellite  Circumnavigation 
%  Main  Program 
%  Capt  Stan  Straight,  USAF 
%  Air  Force  Institute  of  Technology 
%  18  Feb  04 

clc;clear; 

global  thetay  thetaz  gammaO  ro_m  a  Tot_TOF  devr  rdO  epsiO 

%%%%%%%%%%%%%%%%%% 

%  Inputs 

%%%%%%%%%%%%%%%%%% 

%  Rotation  of  the  path  from  the  y-z  plane  about  the  y  axis 

thetayd  =  60;  %  degrees 

thetay  =  thetayd. *pi./l  80;  %  radians 

%  Rotation  of  the  path  about  the  z-axis 

thetazd  =  30;  %  degrees 

thetaz  =  thetazd*pi/180;  %  radians 

%  Define  the  initial  angle  of  the  initial  position  from  the  y-axis 
%  direction 

gammaOd  =  45;  %  degrees 

gammaO  =  gamma0d*pi/180;  %  radians 

%  Input  the  radius  of  the  reference  orbit  in  kilometers 
a  =  6778;  %  km 

%  Input  the  time  of  flight  for  the  complete  circumnavigation  as  a  fraction 
%  of  the  reference  orbit's  period 
Tot_TOF  =  .l; 

%  Define  the  nominal  radius 
ro_m  =  100;  %  m 
rO  =  ro_m./1000;  %  km 

%  Maximum  Deviation  Radius 
devr  =  0.04;  %km 

%  Minimum  Radius  to  search  over  (fraction  of  the  devr) 
minrc  =  1.001; 

%  Establish  the  initial  radius  from  the  nominal  path  and  the  angle  from  the 


Version  6.5,  Release  13.  See  Reference  6. 

^  Only  core  programs  and  sub -programs  shown:  main  program  and  sub-programs  call  other  sub-programs 
for  formatting  data  for  plotting  not  included 
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%  minor  axis 
rdO  =  0; 
epsiO  =  0; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  End  Inputs 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Inner  constraint  radius 
rc  =  rO  -  devr; 
maxrc  =  rc+devr; 

%  Outside  constraint  radius 
rt  =  rO  +  devr; 

%  Calculate  alpha  angle  for  algorithm 
al  =  acos(rc./r0); 

%  Establish  the  search  region 
rd  =  minrc.*rc:. 0001:  maxrc; 

%  Call  the  algorithm  as  a  function  firstorddes 
[states, tot_delvee, const, numbpts]  =  firstorddes  (rt,rc,rd); 

%  Check  for  the  plot  of  the  points  where  there  is  a  change  in  the  number  of 
%  burn  points  between  each  different  case 
combo  =  [rd'  numbpts']; 
jun=  1; 

for  s  =  1  :length(rd)- 1 

if  combo(s,2)  >  combo(s+l,2)  I  combo(s,2)  <  combo(s+l,2); 
ptbptijun,:)  =  combo(s,:); 
jun  =  jun  +1; 
end 
end 

%  Determine  the  minimum  point,  and  the  design  radius  producing  the  minimum 
%  value 

[yuck,ind]  =  min(tot_delvee); 
rdyuck  =  rd(ind); 

%  Add  one  to  the  number  of  burn  points  vector  to  determine  the  acctual 
%  number  of  burn  points  within  a  given  state 
bvec  =  numbpts  +  1 ; 

%  Eind  the  plots  for  the  min  state 
smin  =  states  (:,ind); 
for  g  =  l:length(smin)-3 

if  sum([smin(g)  smin(g+l)])  >  0 
xstmin(g)  =  smin(g); 
end 
end 

%  Place  the  burn  points  for  the  minimum  design  vector 
[bposi,rnom]  =  plCostfcnavl(xstmin); 

%  Use  the  optimization  program  to  evaluate  a  new  minima  if  available 
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[xstop,fvalop]  =  genoptimization(xstmin); 

%  Place  the  burn  points  for  the  optimized  state  vector 
[bposiop,rnomo]  =  plCostfcnavl(xstop); 

%  Create  the  actual  intermediate  points  for  each  state  vector 
cvecmin  =  plotintptsavl(xstmin); 
cvecop  =  plotintptsavl(xstop); 

%  Plot  the  total  impulse  as  a  function  of  the  design  radii 
figure(l) 

plot(rd,tot_del  vee, ' .  ) 

%title( ['Analytical  Approach:  Total  Delta  Vee  vs  Change  in  Design  Radius,  devr=',num2str(devr)]) 
xlabel('r_d  (km)') 
ylabel('\Delta  v_t  (km/s)') 
hold  on 

for  h  =1  :length(ptbpt) 

plot( [ptbpt(h,  1 )  ptbpt(h,  1 )] ,  [min(tot_delvee)  max(tot_delvee)]  ,'r : ') 
end 

hold  off 

%  Plot  the  number  of  burn  points 

figure(2) 

plot(rd,bvec) 

%title('Number  of  Burn  Points  Computed  vs.  Design  Radius') 
xlabel('r_d  (km)') 

ylabel('Number  of  Burn  Points  (Includes  the  Initial  Point)') 

%  Call  a  function  to  create  the  torus  surface  representing  the  maximum 
%  deviation  constraint 

[rdx,rdy,rdz]  =  torusmakerl(rO,devr,thetay,thetaz); 
figure(3) 

plot3(bposi(:,l),bposi(:,2),bposi(:,3),'ro',cvecmin(:,l),cvecmin(:,2),cvecmin(:,3),'r-',... 
bposiop(:,l),bposiop(:,2),bposiop(:,3),'k+',cvecop(:,l),cvecop(:,2),cvecop(:,3),'k--',... 
rnom(:,l),rnom(:,2),rnom(:,3),'-.') 
grid  on 

%  title( ['Optimized  Flight  Path,  Thetay=',num2str(thetayd),',  Thetaz=',num2str(thetazd),... 

%  ',  Gamma0',num2str(gamma0d),',  b=',num2str(numbpts(ind)+l),',  TOF=',num2str(Tot_TOF)]) 

legend('Design  Burn  Points', 'Design  Path', 'Optimized  Burn  Points', 'Optimized  Path', 'Nominal  Path') 

xlabel('x  (km)') 

ylabelCy  (km)') 

zlabel('z  (km)') 

axis  square 

hold  on 

mesh(rdx,rdy,rdz) 
hidden  off 
shading  flat 

colormap([0.9  0.9  0.9]) 
hold  off 
axis  equal 

%  Create  the  plot  showing  the  deviation  from  the  nominal  path  for  each  the 
%  design  state  and  the  optimized  state  vector 
[cmi,tvecmi]  =  planonlincont(xstmin); 
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[co,tveco]  =  planonlincont(xstop); 
cmin  =  zero  s(length(tvecmi)+ 1,1)'; 
tvecmin  =  zeros(length(tvecmi)+l,l)'; 
cmin(2:length(tvecmi)+l)  =  cmi; 
tvecmin(2  :length(tvecmi)+ 1  )=  tvecmi; 

cop  =  zeros(length(tvecmi)+l,l)'; 
tvecop  =  zeros(length(tvecmi)+l,l)'; 
cop(2:length(tvecmi)+l)  =  co; 
tvecop(2 :  length(tvecmi)+ 1  )=  tveco ; 

%  Compute  a  line  for  visualizing  the  constraint 
dep  =  ones(length(tvecmi)+l,l)'; 
devrplot  =  dep.*devr; 

figure(4) 

plot(tvecmin, cmin,'.-’, tvecop, cop, tvecmin, devrplot,’-') 

%  title( ['Deviation  of  Actual  Path','Thetay=',num2str(thetayd),',  Thetaz=',num2str(thetazd),... 

%  ',  Gamma0',num2str(gamma0d),',  b=',num2str(numbpts(ind)-i-l),',  TOF=',num2str(Tot_TOF)]) 

xlabel('Time  (fraction  of  TOF)') 
ylabel('\rho  (km)') 
grid  on 

axis([0  1  0  devr-i-devr.*.05]) 

legend('Design  Path', 'Optimized  Path', 'Constraint  Boundary',4) 
b0=  (length(xstmin)-i-2)./4; 

%  Create  a  dummy  variable  to  populate  with  the  EAET  method 
xst0(l:2.*b0)  =  zeros(2.*b0,l)'; 

%  Equal  Angle/Equal  Time 
gammai  =  1  ./bO; 
ti  =  (l./bO); 
fori=  I:b0-1 

xst0(2.*b0  -I-  i)  =gammai; 
xst0(2.*b0  -I-  i-i-(bO- 1))  =  ti; 
end 

%  Compute  the  total  impulse  for  the  EAET 
eqtot_delvee  =  Costfcnavl(xstO); 

disp(['Equal  Burn/Equal  Time  Total  Delta  Vee:  ',num2str(eqtot_delvee)]) 
disp(['Design  Total  Delta  Vee:',num2str(yuck)]) 
disp(['  Found  at  an  rd  of:',num2str(rdyuck)]) 
disp( ['Optimized  Total  Delta  Vee:',num2str(fvalop)]) 

%  Compute  the  percentage  difference  between  the  Design  and  Optimized 
%  Total  Delta  Vee 

designper  =  ((eqtot_delvee  -  yuck)./eqtot_delvee).*100; 
optimper  =  ((eqtot_delvee  -  fvalop)./eqtot_delvee).*100; 
dispC  ') 

disp( ['Percentage  Saved  by  Design:',num2str(designper),'%']) 
disp( ['Percentage  Saved  by  Optimization:',num2str(optimper),'%']) 
if  const(ind)  >  0 

disp('*******Designed  State  Exceeds  Constraints*******') 
end 
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Subprogram:  Analytical  Design  Function 

function  [states, tot_delvee, const, numbpts]  =  firstorddes(rt,rc,rd) 

%  This  function  is  the  analytical  design  method  given  the  inner 
%  contraint  radius,  rc,  the  last  point  radius  rt,  and  the  search  region  for 
%  the  design  radius  (rd).  The  design  radius  must  be  larger  than  one  point. 
%  The  output  is  the  states  for  each  design  radius,  the  total  impulse  for 
%  those  states  and  the  constraint  matrix  for  those  states.  If  any  number 
%  in  the  constraint  matrix  is  greater  than  0,  then  the  state  exceeds  the 
%  constraints.  The  format  for  calling  this  function  is 
%  [states,tot_delvee,const]  =  firstorddes(rt,rc,rd) 

global  thetay  thetaz  gammaO  ro_m  a  Tot_TOF  devr  rdO  epsiO 

rO  =  ro_m./1000;  %  km 
maxrc  =  rc+devr; 
al  =  acos(rc./r0); 
sigma  =  0; 


for  k  =  l:length(rd) 

%  First  path  length 

p(l,k)  =  sqrt(r0^2-rc^2)+sqrt(rd(k)^2-rc^2); 
pf(k)  =  sqrt(rF^2-rc'^2)+sqrt(rd(k)'^2-rc'^2); 

gamma(l,k)  =  asin((p(l,k).*rc)./(rO.*rd(k))); 
gammaf(k)  =  asin((pf(k).*rc)./(rt.*rd(k))); 

counter  =  1 ; 

while  sigma(k)  <  (2.*pi-gammaf(k)) 

gamma(counter+l,k)  =  2.*acos(rc./rd(k)); 
counter  =  counter  +  1; 
sigmat  =  sum(gamma,l); 
sigma(k)  =  sigmat(k); 
end 

if  sigma(k)  >  (2.*pi  -  gammaf(k)) 

gamt(l:  counter- l,k)  =  gamma(l :  counter- l,k); 
sigmaft  =  sum(gamt,l); 
sigmaf(k)  =  sigmaft(k); 

gamma(counter,k)  =  2.*pi  -  sigmaf(k)  -  gammaf(k); 
end 

sigma(k-Hl)  =  0; 
end 

%  Normalize  the  gamms  to  2  pi 
gamman  =  gamma./(2.*pi); 

%  Set  the  time  between  the  burn  points  to  the  same  ratio  with  respect  to 
%  the  total  time 
timen  =  gamman; 

%  Compute  the  normalized  distance  from  the  nominal 
rdn  =  abs(rO-rd)./devr; 
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%  Set  all  the  points  except  for  the  last  one  on  the  inside  of  the  nominal 
%  within  the  nominal’s  plane;  place  the  last  point  on  the  outside  of  the 
%  torus  in  the  nominal's  plane 
epsi(l:length(rd)-l)  =  0.5; 
epsi(length(rd))  =  0; 

%  Build  each  unique  state  vector 
gut  =  size(gamman); 

for  h  =  l:length(rd) 
for  q  =  l:gut(l) 
if  gamman(q,h)>0 
numbpts(h)  =  q; 
hi  =  numbpts; 
end 
end 
end 

%  Add  one  to  the  number  of  burn  points  vector  to  determine  the  acctual 
%  number  of  burn  points  within  a  given  state 
bvec  =  numbpts  +  1 ; 
for  u  =  1  :length(rd) 
xstn(l:bi(u))  =  rdn(u); 
xstn(bi(u)+l)  =  abs(rO-rt)./devr; 
xstn(bi(u)+2:2.*bi(u)+l)  =  0.5; 
if  rt  >=  rO 

xstn(2.*bi(u)+2)  =  0; 
else 

xstn(2.  *bi(u)+2)=0.5 ; 
end 

xstn(2.*bi(u)+3:3.*bi(u)+2)  =  gamman(l :  numbpts  (u),u); 
xstn(3.*bi(u)+3:4.*bi(u)+2)  =  timen(l:numbpts(u),u); 
if  u  ==  1 

states  =  zeros(length(xstn),length(rd)); 
end 

xstn  =  xstn'; 

states(l:length(xstn),u)  =  xstn; 

%  Compute  the  cost  function  value  for  the  new  state  vector 
tot_delvee(u)  =  Costfcnavl(xstn); 

%  Compute  the  constraint  function  for  the  state  to  determine  if  the 
%  constraints  are  violated 
[dum,ceq]  =  anonlincon(xstn); 
c(u,l:length(dum))  =  dum; 
const(u)=0; 
constc(u)  =  0; 
for  y  =  1  :length(c) 
if  c(u,y)  >  0 
const(u)  =  const(u)+l; 
constc(u)  =  1 ; 
end 
end 

xstn  =  0; 
dum  =  []; 
end 
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Sub-Subprogram:  Cost  Function  Evaluation  -  Analytical  Design 


function  tot_delvee  =  costfcnavl(xst) 

%  This  is  the  cost  function  to  be  used  with  fmincon.  It  computes  the  total 
%  delta  vee  required  to  perform  a  circumnavigation,  global  variables  are  used  and 
%  should  be  defined  by  the  main  program,  additionally  the  state  vector  xst  should  be 
%  of  the  format  rd's,  epsilon's,  gamma's  then  time  and  of  4*b  -2  length. 

global  thetay  thetaz  gammaO  ro_m  a  Tot_TOF  devr  rdO  epsiO 

%  Calculate  the  mean  motion  of  the  reference  orbit  in  rad/s 
n  =  sqrt((3.98601*10^5)/a^3); 

%  Compute  the  time  of  circumnavigation 
refperiod  =  2*pi/n;  %  seconds 
tot_toc  =  Tot_TOF*refperiod;  %  seconds 

b=  (length(xst)+2)./4; 
bmax  =  b; 

%  Path  radius  in  kilometers 
rO  =  ro_m/1000;  %  kilometers 

%  Calculate  the  real  values  from  the  normalized  states 

rds  =  xst(l:b).*devr; 

epsis  =  xst(b+l:2.*b).*2.*pi; 

gamms  =  xst(2.*b+l:3.*b-l).*2.*pi; 

tims  =  xst(3.*b:4.*b-2).*tot_toc; 

%  Define  the  initial  position  vector  of  the  initial  burn  point 
rdxO  =  rdO.*sin(epsiO).*cos(thetay).*cos(thetaz)  + 

(rO+rdO.  *cos(epsiO)) .  *  sin(gammaO) .  *cos(thetaz) .  *  sin(thetay) . . . 

-  (rO+rdO.  *cos(epsiO)) .  *cos(gammaO) .  *  sin(thetaz) ; 
rdyO  =  (rO+rdO.  *cos(epsiO)) .  *cos(gammaO) .  *cos(thetaz)+rdO.  *  sin(epsD).  *cos(thetay) .  *sin(thetaz). . . 

+(rO+rdO.  *cos(epsiO)) .  *  sin(gammaO) .  *  sin(thetay) .  *  sin(thetaz) ; 
rdzO  =  (rO+rdO.*cos(epsiO)).*sin(gammaO).*cos(thetay)-rdO.*sin(epsiO).*sin(thetay); 

rint  =  [rdxO;rdyO;rdzO]'; 

%  Define  the  initial  velocity  vector  of  the  initial  burn  point 
Vint  =[000]; 

%  Establish  the  initial  velocity  of  the  interceptor  at  the  initial  point 
dvf(l,:)  =  vint; 

for  s  =  1  :b 
if  s  <  b 

rd  =  rds(s); 
epsi  =  epsis(s); 
toe  =  tims(s); 
if  s  ==1 

gamma^j  =  gamms  (s)+gamma0; 
else 

gamma  J  =  gamms(s)+gammaj; 
end 
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elseif  s  =  b 
rd  =  rds(s); 
epsi  =  epsis(s); 

gammaj  =  (2.*pi-  suni(gamms))+gammaj; 
toe  =  tot_toc  -  sum(tims); 
end 

%Position  of  the  next  burn  point 

rdx  =  rd.*sin(epsi).*cos(thetay).*cos(thetaz)  + 

(rO+rd.  *cos(epsi)) .  *  sin(gamma_j) .  *cos(thetaz) .  *  sin(thetay) . . . 
-(rO+rd.*cos(epsi)).*cos(gamma_j).*sin(thetaz); 
rdy  =  (rO+rd.  *cos(epsi)) .  *cos(gamma_j ) .  *cos(thetaz)+rd.  *  sin(epsi) .  *cos(thetay) .  *  sin(thetaz) . . . 

+(rO+rd.  *cos(epsi)) .  *  sin(gamma_j ) .  *  sin(thetay) .  *  sin(thetaz) ; 
rdz  =  (rO+rd.  *cos  (epsi)) .  *  sin(ganinia  J ) .  *  co  s  (thetay  )-rd.  *  sin(epsi) .  *  sin(thetay) ; 

rfin  =  [rdx;rdy;rdz]'; 

%  Compute  the  velocity  needed  to  go  from  the  current  point  to  the  next  point 
%  in  the  given  time  of  flight 
dvi(s,:)  =  hills  vel2(rint, rfin, toe, n); 
if  s  <  b 

%  Compute  the  actual  velocity  at  the  next  point 
dvf(s+l,:)  =  hillsvelf(rint,dvi(s,:),toc,n); 

end 

%  Reset  the  position  of  the  next  burn  point  to  the  current  burn 
%  point  to  propagate  the  next  sequential  burn  point  along  the  path 
rint  =  rfin; 


end 

%  Compute  the  change  in  velocity 
delvee_vec  =  dvf  -  dvi; 

if  b  ==  1 

delvee_mag_burn(b)  =  norm(delvee_vec(b,:)); 
end 

for  g  =  1  :b 

delvee_mag_burn(g)  =  norm(delvee_vec(g,:)); 
end 

tot_delvee  =  sum(delvee_mag_burn); 


Sub-Subprogram:  Compute  Impulse  Between  Points  Using  Hill's  Equations 

function  vint  =  hillsvel2(rint,rfin,tof,n) 

%  This  program  will  take  an  initial  and  final  position  row  vectors  defined 
%  in  the  LVLH  frame,  the  time  of  flight  between  the  points  and  output  the 
%  initial  velocity  required  to  achieve  these  parameters.  The  output  is  in 
%  a  row  vector  format  for  the  velocity  components  in  the  LVLH  frame.  The 
%  time  of  flight  needs  to  be  in  seconds  and  the  value  of  n  has  to  be  in 
%  rad/s.  n  is  the  mean  motion  of  the  reference  orbit. 
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%  Define  the  value  of  psi 
ang  =  n  *  tof; 

%  Input  the  phi  rr  and  phi  rv  matrices 

phirr  =  [(4-3*cos(ang)),0,0;6*(sin(ang)-ang),l,0;0,0,cos(ang)]; 

inphirv  =  [n*(3*ang-4*sin(ang))/(-8+8*cos(ang)+3*ang*sin(ang)),-2*n*(cos(ang)-l)/(- 
8+8*cos(ang)+3*ang*sin(ang)),0; 

2*n*(cos(ang)-l)/(-8+8*cos(ang)+3*ang*sin(ang)),n*sin(ang)/(8-8*cos(ang)-3*ang*sin(ang)),0; 

0,0,n/sin(ang)]; 

%  compute  the  required  initial  velocity 
vint  =  (inphirv  *(rfin'-phirr*rint'))'; 

%end 


Sub-Subprogram:  Propagate  Velocity  from  Previous  Burn  Point  with  HiWs 
Equations 

function  vfin  =  hillsvelf(rint, vint, tof, n) 

%  This  program  will  take  an  initial  position  and  initial  velocity  row  vectors  defined 
%  in  the  LVLH  frame,  the  time  of  flight  between  the  points  and  output  the 
%  velocity  at  the  end  of  the  time  of  flight  required  to  achieve  these  parameters. 

%  The  output  is  in 

%  a  row  vector  format  for  the  velocity  components  in  the  LVLH  frame.  The 
%  time  of  flight  needs  to  be  in  seconds  and  the  value  of  n  has  to  be  in 
%  rad/s.  n  is  the  mean  motion  of  the  reference  orbit. 

%  Define  the  value  of  psi 
ang  =  n  *  tof; 

%  Input  the  phi  vr  and  phi  vv  matrices 

phivr  =  [3*n*sin(ang),0,0;6*n*(cos(ang)-l),0,0;0,0,-n*sin(ang)]; 

phivv  =  [cos(ang),2*sin(ang),0;-2*sin(ang),-3+4*cos(ang),0;0,0,cos(ang)]; 

%  compute  the  required  final  velocity 
vfin  =  (phi  vr*rint'+phivv*  vint')'; 

Subprogram:  Perform  Optimization  Using /mincon 


%  Compute  the  optimized  path  given  the  global  constraints  and  initial  guess 
%  state  vector 
%  Capt  Stan  Straight 

function  [xst,fval]  =  genoptimization(xstO) 

global  thetay  thetaz  gammaO  ro_m  a  Tot_TOF  devr  rdO  epsiO 
%  Establish  the  options  parameters  for  the  optimization  function 
dmaxchange  =  le-1; 
dminchange  =  le-9; 
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TolFun  =  2e  -9; 

TolCon  =  le-9; 

TolX  =  le-9; 

%  Calculate  the  mean  motion  of  the  reference  orbit  in  rad/s 
n  =  sqrt((3.98601*10^5)./a^3); 

%  Compute  the  time  of  circumnavigation 
refperiod  =  2*pi./n;  %  seconds 
tot_toc  =  Tot_TOF*refperiod;  %  seconds 
%  Compute  the  number  of  burn  points 
b=  (length(xst0)-i-2)./4; 
bmax  =  b; 

%  Path  radius  in  kilometers 
rO  =  ro_m/1000;  %  kilometers 

%  Establish  the  A  matrix  contraint  and  the  b  matrix  constraint 
A  =  zeros(4.*b-2,4.*b-2); 

A(l,2.*b-rl:3.*b-l)  =  ones(l,b-l); 

A(2,3.*b:4.*b-2)  =  ones(l,b-l); 

%  Establish  the  b  matrix  contraint 
bineq  =  zeros(4.*b-2,l); 
bineq(l)  =  1; 
bineq(2)  =  1 ; 

%  Define  the  lower  and  upper  bounds  of  the  functions 
for  k  =  l:4.*b-2 
ub(k)  =  1; 
end 

ub=ub’; 

for  p  =  l:4.*b-2 
lb(p)  =  0; 
end 

lb(3.*b:4.*b-2)  =  le-7; 
lb  =  lb’; 
options  = 

optimset('LargeScale', 'off, 'MaxFunEvals', 20000, 'Maxiter', 500, 'TolFun', TolFun, 'TolCon', TolCon,... 
'TolX',TolX, 'Display', 'iter', 'DiffMaxChange',dmaxchange,'DiffMinChange',dminchange); 

%  Use  fmincon  function  to  find  the  state  vector  that  produces  the  minimum  total  deltavee 
[xst,fval,exitflag,output,lambda,grad,hessian]  = 

fmincon(  @  Costfcngv  1  ,xst0,  A,bineq,  [] ,  []  ,lb,ub,  @  gnonlincon,options) ; 


Sub-Subprogram:  Cost  Function  Evaluation  for  Optimization. 

function  tot_delvee  =  costfcngv  l(xst) 

%  This  is  the  cost  function  to  be  used  with  fmincon.  It  computes  the  total 
%  delta  vee  required  to  perform  a  circumnavigation,  global  variables  are  used  and 
%  should  be  defined  by  the  main  program,  additionally  the  state  vector  xst  should  be 
%  of  the  format  rd's,  epsilon's,  gamma's  then  time  and  of  4*b  -2  length. 

global  thetay  thetaz  gammaO  ro_m  a  Tot_TOF  devr  rdO  epsiO 

b=  (length(xst)-i-2)./4; 
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bmax  =  b; 


%  Calculate  the  mean  motion  of  the  reference  orbit  in  rad/s 
n  =  sqrt((3.98601*10^5)/a^3); 

%  Compute  the  time  of  circumnavigation 
refperiod  =  2*pi/n;  %  seconds 
tot_toc  =  Tot_TOF*refperiod;  %  seconds 

%  Path  radius  in  kilometers 
rO  =  ro_m/1000;  %  kilometers 

%  Calculate  the  real  values  from  the  normalized  states 

rds  =  xst(l:b).*devr; 

epsis  =  xst(b+l:2.*b).*2.*pi; 

gamms  =  xst(2.*b+l:3.*b-l).*2.*pi; 

tims  =  xst(3.*b:4.*b-2).*tot_toc; 

%  Define  the  initial  position  vector  of  the  initial  burn  point 
rdxO  =  rdO.*sin(epsiO).*cos(thetay).*cos(thetaz)  + 

(rO+rdO.  *cos(epsiO)) .  *  sin(gammaO) .  *cos(thetaz) .  *  sin(thetay) . . . 

-  (rO+rdO.  *cos(epsiO)) .  *cos(gammaO) .  *  sin(thetaz) ; 

rdyO  =  (rO+rdO.*cos(epsiO)).*cos(gammaO).*cos(thetaz)+rdO.*sin(epsiO).*cos(thetay).*sin(thetaz)... 

+(rO+rdO.*cos(epsiO)).*sin(gammaO).*sin(thetay).*sin(thetaz); 
rdzO  =  (rO+rdO.*cos(epsiO)).*sin(gammaO).*cos(thetay)-rdO.*sin(epsiO).*sin(thetay); 

rint  =  [rdxO;rdyO;rdzO]'; 

%  Define  the  initial  velocity  vector  of  the  initial  burn  point 
vint  =[000]; 

%  Establish  the  initial  velocity  of  the  interceptor  at  the  initial  point 
dvf(l,:)  =  vint; 

for  s  =  1  :b 
if  s  <  b 
rd  =  rds(s); 
epsi  =  epsis(s); 
toe  =  tims(s); 
if  s  =1 

gammaj  =  gamms(s)+gammaO; 
else 

gammaj  =  gamms(s)+gammaj; 
end 

elseif  s  ==  b 
rd  =  rds(s); 
epsi  =  epsis(s); 

gammaj  =  (2.*pi-  sum(gamms))+gammaj; 
toe  =  tot_toc  -  sum(tims); 
end 

%Position  of  the  next  burn  point 

rdx  =  rd.*sin(epsi).*cos(thetay).*cos(thetaz)  + 

(rO+rd.  *cos(epsi)) .  *  sin(gamma  J ) .  *cos(thetaz) .  *  sin(thetay) . . . 

-  (rO+rd.  *cos(epsi))  .*  cos  (gamma J) .  *  sin(thetaz) ; 
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rdy  =  (rO+rd.  *cos(epsi)) .  *cos(gamma_j) .  *cos(thetaz)+rd.  *sin(epsi).  *cos(thetay) .  *sin(thetaz) . . . 

+(rO+rd.  *cos(epsi)) .  *  sin(gamma  J ) .  *  sin(thetay) .  *  sin(thetaz) ; 
rdz  =  (rO+rd.  *cos(epsi)) .  *  sin(gamma  J) .  *cos(thetay)-rd.  *sin(epsi).  *sin(thetay) ; 

rfin  =  [rdx;rdy;rdz]'; 

%  Compute  the  velocity  needed  to  go  from  the  current  point  to  the  next  point 
%  in  the  given  time  of  flight 
dvi(s,:)  =  hills  vel2(rint, rfin, toe, n); 
if  s  <  b 

%  Compute  the  actual  velocity  at  the  next  point 
dvf(s+l,:)  =  hillsvelf(rint,dvi(s,:),toc,n); 

end 

%  Reset  the  position  of  the  next  burn  point  to  the  current  burn 
%  point  to  propagate  the  next  sequential  burn  point  along  the  path 
rint  =  rfin; 

end 

%  Compute  the  change  in  velocity 
delvee_vec  =  dvf  -  dvi; 

if  b  ==  1 

delvee_mag_burn(b)  =  norm(delvee_vec(b,:)); 
end 

for  g  =  1  :b 

delvee_mag_burn(g)  =  norm(delvee_vec(g,:)); 
end 

tot_delvee  =  sum(delvee_mag_burn); 


Sub-Subprogram:  Determine  Flight  Path  Constraints 

function  [c,ceq]  =  gnonlincon(xst) 

global  thetay  thetaz  gammaO  ro_m  a  Tot_TOF  devr  rdO  epsiO 

%  Define  the  number  of  intermediate  points 
z  =  20; 

%  Calculate  the  mean  motion  of  the  reference  orbit  in  rad/s 
n  =  sqrt((3.98601*10^5)/a^3); 

%  Compute  the  time  of  circumnavigation 
refperiod  =  2*pi/n;  %  seconds 
tot_toc  =  Tot_TOF*refperiod;  %  seconds 

%  Path  radius  in  kilometers 
rO  =  ro_m/1000;  %  kilometers 
b=  (length(xst)+2)./4; 
bmax  =  b; 


%  Calculate  the  real  values  from  the  normalized  states 

rds  =  xst(l:b).*devr; 

epsis  =  xst(b+l:2.*b).*2.*pi; 
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gamms  =  xst(2.*b+l:3.*b-l).*2.*pi; 
tims  =  xst(3.*b:4.*b-2).*tot_toc; 

%  Define  the  initial  position  vector  of  the  initial  burn  point 
rdxO  =  rdO.*sin(epsiO).*cos(thetay).*cos(thetaz)  + 

(rO+rdO.  *cos(epsiO)) .  *  sin(ganiniaO) .  *cos(thetaz) .  *  sin(thetay) . . . 

-  (rO+rdO.  *cos(epsiO)) .  *cos(ganiniaO) .  *  sin(thetaz) ; 

rdyO  =  (rO+rdO.*cos(epsiO)).*cos(gammaO).*cos(thetaz)+rdO.*sin(epsiO).*cos(thetay).*sin(thetaz)... 

+(rO+rdO.  *cos(epsiO)) .  *  sin(gammaO) .  *  sin  (thetay) .  *  sin(thetaz) ; 
rdzO  =  (rO+rdO.*cos(epsiO)).*sin(gammaO).*cos(thetay)-rdO.*sin(epsiO).*sin(thetay); 

rint  =  [rdxO;rdyO;rdzO]’; 

%  Define  the  initial  velocity  vector  of  the  initial  burn  point 
vint  =[000]; 

%  Establish  the  initial  velocity  of  the  interceptor  at  the  initial  point 
dvf(l,:)  =  vint; 

for  s  =  1  :b 
if  s  <  b 

rd  =  rds(s); 
epsi  =  epsis(s); 
toe  =  tims(s); 
if  s  =1 

gammaj  =  ganinis(s)+ganiniaO; 
else 

gammaj  =  gamms(s)+gammaj; 
end 

elseif  s  =  b 
rd  =  rds(s); 
epsi  =  epsis(s); 

gammaj  =  (2.*pi-  sum(gamms))+gamma J ; 
toe  =  tot_toc  -  sum(tims); 
end 

%Position  of  the  next  burn  point 

rdx  =  rd.*sin(epsi).*cos(thetay).*cos(thetaz)  + 

(rO+rd.  *cos(epsi)) .  *  sin(gamma  J) .  *cos(thetaz) .  *  sin(thetay) . . . 

-  (rO+rd.  *cos(epsi))  .*  cos  (gamma J) .  *  sin(thetaz) ; 

rdy  =  (rO+rd.*cos(epsi)).*cos(gammaJ).*cos(thetaz)+rd.*sin(epsi).*cos(thetay).*sin(thetaz)... 

+(r0+rd.  *cos(epsi)) .  *  sin(gamma  J ) .  *  sin(thetay) .  *  sin(thetaz) ; 
rdz  =  (rO+rd.  *cos(epsi)) .  *  sin(gamma  J) .  *cos(thetay)-rd.  *sin(epsi).  *sin(thetay) ; 

rfin  =  [rdx;rdy;rdz]'; 

%  Compute  the  velocity  needed  to  go  from  the  current  point  to  the  next  point 
%  in  the  given  time  of  flight 
dvi(s,:)  =  hills  vel2(rint, rfin, toe, n); 

%  Compute  the  vector  of  intermediate  points  within  the 
cirfin  =  inthillsoptnonlin(rint,dvi(s,:),toc,n,z); 
q(s)  =  (z+l+(s-2).*z)’; 
cvec(q(s):q(s)+z-l,:)  =  cirfin; 
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%  Reset  the  position  of  the  next  burn  point  to  the  current  burn 
%  point  to  propagate  the  next  sequential  burn  point  along  the  path 
rint  =  rfin; 


end 

%  Create  a  rotation  matrix  from  the  LVLH  frame  to  the  circumnav  frame  given 
%  the  values  of  gammaO,  thetay,  and  thetaz 

rotxtop  =  [cos(thetay) .  *cos(thetaz)  ,cos(thetay) .  *  sin(thetaz),  - 1 .  *  sin(thetay) ; . . . 
cos(thetaz) .  *  sin(gammaO) .  *  sin(thetay  )-cos(gammaO) .  *  sin(thetaz) , . . . 

cos(gammaO) .  *cos(thetaz)+sin(thetay) .  *  sin(thetaz) .  *  sin(gammaO),cos(thetay) .  *  sin(gammaO) 
cos(gammaO) .  *cos(thetaz) .  *  sin(thetay)+sin(gammaO) .  *  sin(thetaz), . . . 

- 1 .  *  (cos (thetaz) .  *  sin(gammaO))+cos  (gammaO) .  *  sin(thetay) .  *  sin(thetaz) , . . . 
cos(gammaO).*cos(thetay)]; 

pathvec  =  zeros(length(cvec),3); 

%  Create  the  c  matrix  by  determining  the  distance  of  the  intermediate 
%  position  to  the  position  of  the  nominal  path, 
for  m  =  l:length(cvec) 

%  Extract  each  intermediate  position 
cvecp  =  cvec(m,:)'; 

%  Rotate  each  vector  into  the  pqw,  then  take  the  projection  onto  the 
%  pq  frame 

cvecprot(m,:)  =  (rotxtop*  cvecp)'; 
cvecproj(m,:)  =  cvecprot(m,2:3); 

pathvec(m,2:3)  =  rO.*(cvecproj(m,:)./norm(cvecproj(m,:))); 
end 

cmat  =  cvecprot-pathvec; 

for  w  =  l:length(cvec) 

c(w)  =  norm(cmat(w,:))-devr; 
end 

%  By  setting  this  quantity  to  zero  it  eliminates  the  nonlinear  equality  constraint 
%  requirement 
ceq  =  0; 
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